Abstract
High-grade serous cancer (HGSC) is the most common subtype of ovarian cancer. HGSC is highly aggressive with poor patient outcomes, and a deeper understanding of HGSC tumorigenesis could help guide future treatment development. To systematically characterize the underlying pathologic mechanisms and intratumoral heterogeneity in human HGSC, we used an optimized single-cell multiomics sequencing technology to simultaneously analyze somatic copy-number alterations (SCNA), DNA methylation, chromatin accessibility, and transcriptome in individual cancer cells. Genes associated with interferon signaling, metallothioneins, and metabolism were commonly upregulated in ovarian cancer cells. Integrated multiomics analyses revealed that upregulation of interferon signaling and metallothioneins was influenced by both demethylation of their promoters and hypomethylation of satellites and LINE1, and potential key transcription factors regulating glycolysis using chromatin accessibility data were uncovered. In addition, gene expression and DNA methylation displayed similar patterns in matched primary and abdominal metastatic tumor cells of the same genetic lineage, suggesting that metastatic cells potentially preexist in the subclones of primary tumors. Finally, the lineages of cancer cells with higher residual DNA methylation levels and upregulated expression of CCN1 and HSP90AA1 presented greater metastatic potential. This study characterizes the critical genetic, epigenetic, and transcriptomic features and their mutual regulatory relationships in ovarian cancer, providing valuable resources for identifying new molecular mechanisms and potential therapeutic targets for HGSC.
Integrated analysis of multiomic changes and epigenetic regulation in high-grade serous ovarian cancer provides insights into the molecular characteristics of this disease, which could help improve diagnosis and treatment.
Introduction
Ovarian cancer is one of the leading causes of gynecologic cancer death (1), due in part to its high heterogeneity. High-grade serous cancer (HGSC) is the most aggressive and common subtype of ovarian cancer (comprising approximately 70% of all cases). HGSC is generally diagnosed at an advanced stage, presenting with extensive metastases and massive ascites. Despite an initial clinical response to the surgery and chemotherapy, the majority of patients will relapse with fatal outcomes. Hence, understanding the pathologic mechanisms of HGSC initiation, promotion, and progression will be of great value.
Intratumor heterogeneity (ITH) is one of the most challenging obstacles to study HGSC (2). Previous studies have revealed the heterogeneity of HGSC in single-nucleotide variations (SNV), somatic copy-number alterations (SCNA), aberrant transcriptomic programs, and tumor subtypes (3–5). However, these studies are usually based on bulk sequencing, which permits only tissue-level resolution and has limited ability to illustrate ITH. In addition, tumor development involves not only alterations in the genome and transcriptome, but also extensive epigenetic changes, which remain to be comprehensively analyzed in HGSC.
Single-cell sequencing provides a powerful tool for the study of ITH. Most studies have used single-cell RNA sequencing (scRNA-seq) to characterize different phenotypic cell states of cancer cells (6–9). However, to accurately explore the dynamic changes of tumor promotion (i.e., the changes between normal tissues and primary tumors) and progression (i.e., the changes between primary tumors and metastases), reconstruction of the genetic lineages of cancer cells is necessary. Naturally occurring genomic mutations, like SNVs and SCNAs, leave the real signatures to trace evolutionary histories of cancer cells, which can be detected by single-cell DNA sequencing. Besides tracing the genetic lineages and revealing the phenotypic changes during tumorigenesis, we also wondered the epigenetic alterations and their regulation on phenotypes during tumor promotion and progression. Thus, we improved single-cell multiomics sequencing technology (scCOOL-seq; ref. 10) to simultaneously assess the SCNAs, DNA methylome, chromatin accessibility, and transcriptome in the same individual cells. Using improved scCOOL-seq, we reconstructed the genetic lineages using SCNAs. Tracing the genetic lineages, we characterized the multiomics changes as well as the epigenetic regulation of phenotypes from normal cells to primary tumor cells and further to matched metastatic cells. Our study offers deeper insights for HGSC and provides a valuable resource for developing new therapeutic strategies.
Materials and Methods
Study approval
This study was approved by the ethics committee of the Peking University People's Hospital (License No. 2019HPB034-01) and was conducted in accordance with the Declaration of Helsinki. All patients provided the written informed consent before the surgery.
Clinical samples
No patient received any tumor-related treatment before surgery. Fallopian tubes (FT) were obtained from patients who underwent hysterectomy with bilateral salpingo-oophorectomy (BSO) for benign and malignant gynecologic diseases (patient OC20, OC21, OC25, OC26, and OC28) and risk-reducing salpingo-oophorectomy (RRSO) surgery (patient OC22). Only epithelial cells obtained from pathologically normal FT tissues and proved to have no SCNAs using improved scCOOL-seq were considered as normal control cells.
Single cell preparation
The fresh specimens were digested into single-cell suspensions using the MACS Tissue Dissociation Kit (Miltenyi Biotec, #130-095-929). Epithelial cells of tumor tissues were isolated by magnetic activated cell sorting (MACS; CD326 Microbeads, Miltenyi Biotec, #130-061-101, RRID: AB_2832928). The live normal fallopian tube epithelial (FTE) cells were sorted by FACS using the anti-human CD326 antibody (BioLegend, #324208, RRID: AB_756082) and 7AAD (BD Pharmingen, #559925).
Library construction of improved scCOOL-seq
The sorted single cells were first treated by the GpC methyltransferase (NEB, #M0227L). Next, each cell was incubated with prewashed carboxylic acid dynabeads (Invitrogen, #65011) to separate the nuclei from the RNA-containing cytoplasm according to the scTrio-seq2 protocol (11).
The supernatant containing mRNA was processed following the modified STRT-seq protocol (12). The nuclear DNA fraction was processed following a previously described method (13). Briefly, the nuclei with magnetic Dynabeads were treated with bisulfite (Zymo Research, #D5044) and were then amplified using random primers with anchor sequence. After purification by AMPure XP beads, the products were amplified by index primers for 16 cycles. The libraries were finally sequenced on the Illumina HiSeq 4000 platform.
Gene knockout, gene overexpression, and cell viability assay
The human ovarian cancer cell line SKOV3 (RRID: CVCL_0532) and ES2 (RRID: CVCL_3509) were purchased from ATCC. Cell lines were authenticated by short tandem repeat (STR) analysis and were mycoplasma negative. Cells within 20 passages were used in the study. We used the clustered regularly interspaced short palindromic repeats (CRISPR)-Cas9 to knockout CCN1. The guide RNA (gRNA) sequences are as follows: gRNA1: CGCGCACTTGGGCGCCTCCA; gRNA2: GCAGATCCCCTTCAGAGCGG. Overexpression of HSPA6 was performed by transfecting the overexpression plasmid using Lipofectamine 3000 Reagent (Invitrogen, #L3000150) as described previously (14). The HSP90AA1 inhibitor TAS-116 was purchased from Selleck Chemicals (#S7716). Cell viability was analyzed by the CellTiter-Glo 2.0 Assay (Promega, #G9241).
Wound healing assay
When cells grew to full confluence, a wound was scratched using a 200-μL micropipette tip. The cells were washed with PBS to remove the cell debris and then grown in serum-free medium. The wound healing results were recorded using a phase-contrast microscope (Leica), and the migration areas were calculated using the ImageJ software.
Transwell assays
Transwell assays were performed with transwell chambers (Coring, #3422) with or without coating with Matrigel (Corning, #356231) to measure the invasion and migration ability, respectively. Cells were harvested in serum-free medium and seeded into the upper chambers. After 24-hour incubation, the cells in the lower chambers were fixed and stained, and cells in the upper chambers were removed with a cotton swab. The photographs were taken using a microscope and the cell numbers were counted using ImageJ software.
Data availability
The raw data have been deposited in the Genome Sequence Archive (GSA) under the accession number HRA000360 and are available upon reasonable request. The processed data associated with SCNAs, DNA methylome, chromatin accessibility, and transcriptome have been submitted to the Gene Expression Omnibus (GEO) database under the accession number GSE189955.
Code availability
Customized codes for data analysis are available at https://github.com/hlxie.
Results
Single-cell multiomics sequencing reveals tumor heterogeneity in HGSC
Here, using improved scCOOL-seq, we profiled normal fallopian tubes, primary tumors, and matched metastases to unveil the molecular characteristics associated with tumor development (Fig. 1A; Supplementary Fig. S1A). To cover ITH as much as possible, the tumor tissues were processed with a multiregional sampling strategy. To avoid variations resulting from treatment effects and disease stages, only patients who were treatment-naïve and had stage III or IV disease were enrolled in our study (Supplementary Table S1).
We first sequenced the RNA compartment of single cells and obtained single-cell transcriptomes for a total of 5,329 cells (Supplementary Table S1). On average, 4,701 genes were detected per cell (Supplementary Fig. S1B). We distinguished cancer cells and normal FTE cells from stromal cells in three steps. First, according to clustering analysis and well-known cell type-specific marker gene expression, cells can be classified as epithelial cells, T cells, macrophages, and fibroblasts (Fig. 1B; Supplementary Fig. S1C). Second, we used the R package InferCNV (15) to infer SCNAs based on averaged expression levels across chromosomal intervals (Fig. 1C). Next, we sequenced the nuclear fractions for some of the epithelial cells from tumor tissues with inferred SCNAs and epithelial cells from FTs without inferred SCNAs. We profiled the SCNA patterns using the DNA methylome data of these selected cells (called “genomic SCNAs” for short; Fig. 1D; Supplementary Figs. S1D and S1E; Supplementary Fig. S2). Finally, we defined the epithelial cells with specific genomic SCNAs as cancer cells and epithelial FT cells without genomic SCNAs as normal FTE cells. Cancer cells display great interpatient heterogeneities as previously reported (6–9), whereas FTE cells and stromal cells exhibit low interpatient heterogeneities (Fig. 1B).
Reconstruction of the genetic lineages of cancer cells
We reconstructed genetic lineages of cancer cells in the same patient using the improved scCOOL-seq genome data. Diverse genetic alterations in cancer cells provide hallmarks for lineage tracing, such as SNVs and SCNAs. Considering that HGSC tumors frequently carry extensive SCNAs, we used SCNAs to construct the genetic lineages. SCNAs were classified into regional aberrations and focal aberrations (16). Focal aberrations frequently occur within chromosome arms and are generally considered irreversible (16, 17). Therefore, we used the breakpoints of subchromosomal focal aberrations to reconstruct the genetic lineages of cancer cells within the same patient (Fig. 1D; Supplementary Fig. S2; ref. 11). Cells containing more than three identical intra-chromosomal breakpoints were considered to belong to the same genetic lineage (11). Our data show that though HGSC carry extensive SCNAs, the composition of genetic lineages in HGSC is relatively simple, compared with colorectal cancer (11). The majority of HGSC patients (8/9) have only one or two lineages of cancer cells (Fig. 1D; Supplementary Fig. S2). Especially for OC16, more than 200 single cells collected from 10 sampling regions all belonged to the same lineage (Supplementary Fig. S2; Supplementary Table S1).
Dynamic changes of RNA expression from normal FTE cells to primary tumor cells
In recent years, the notion that HGSC originates from fallopian tube epithelium has been supported by increasing amounts of evidences (18). Thus, we used FTE cells as controls for the following analyses. According to the genetic lineages, all cancer cells in a given patient originated from the same ancestor cancer cell (Fig. 1D; Supplementary Fig. S2). Thus, we compared all primary tumor cells from the same patient with normal FTE cells to explore the phenotypic changes associated with tumor promotion.
As shown in Fig. 2A, 187 genes were downregulated in primary tumors. Gene ontology (GO) analysis shows that the pathways of eukaryotic translation process and protein ubiquitin were enriched (Fig. 2B), implying great alterations in protein levels. As an important organ to produce and response to hormones, the expression of genes related to hormone responses was also aberrant in cancer cells.
On the other hand, 361 genes were upregulated in primary tumors (Fig. 2A). GO analysis revealed that these differentially expressed genes (DEG) were mainly enriched in pathways associated with VEGFA-VEGFR2 signaling (such as VEGFA, PGK1, and NRP2), metal ion homeostasis (such as MT1E, MT1F, MT1G, MT1H, MT1X, and MT2A), and regulation of viral process (such as ISG15, BST2, STAT1, and IFI27; Fig. 2B). Of note, genes associated with both of oxidative phosphorylation (OXPHOS; such as ATP5MC3, ATP5ME, and ATP6AP1) and aerobic glycolysis (such as ALDOA, ENO1, and LDHA) were upregulated. To investigate whether this phenomenon is due to the existence of two cell subtypes (one upregulated OXPHOS, and the other one upregulated glycolysis), we calculated OXPHOS and glycolysis scores for each cell (see Supplementary Materials and Methods). The FTE cells and cancer cells could be clearly distinguished, whereas cancer cells could not be further separated into subtypes and there were minute variations between patients, suggesting that cancer cells simultaneously elevated the expression of both OXPHOS and aerobic glycolytic genes compared with FTE cells (Fig. 2C–E). Similarly, the patients of TCGA also cannot be divided into two groups according to the expression of OXPHOS and glycolysis (Fig. 2F).
To further clarify the typical characteristics of HGSC, we also compared HGSC with four non-HGSC tumors. We found that HGSC tumors could be clearly separated from non-HGSC tumors by principal component analysis (PCA), demonstrating that HGSC had very distinct transcriptional patterns (Supplementary Figs. S3A and S3B). Furthermore, the differential expression analysis shows that expression of metallothionein genes and interferon signaling pathway genes were also increased in HGSC when compared with non-HGSC tumors (Supplementary Figs. S3C and S3D). Moreover, we performed IHC staining of ISG15 and BST2 to verify this result at protein levels (Fig. 2G; Supplementary Figs. S3E and S3F).
The regulation of SCNAs and DNA methylation on gene expression
HGSC frequently carried strong amplification on chromosome 8 (Chr8; ref. 4). We wondered whether the 361 DEGs upregulated in primary tumors may be resulted from the amplification of Chr8. We identified that 46 genes of these 361 genes were on Chr8, which is the most significantly enriched chromosome for these DEGs (P = 2.7E−14, hypergeometric test; Fig. 3A). Interestingly, these 46 genes were most enriched in the pathways associated with electron transport chain and OXPHOS (Fig. 3B). We also found that the upregulation of these 46 genes was associated with shorter overall survival time (Fig. 3C). In addition, we also observed that MYC (located on Chr8), which regulated cell division, OXPHOS and aerobic glycolysis, was upregulated in 10 out of the 12 patients with HGSC we analyzed (Supplementary Fig. S4A). The results suggested that the amplification of Chr8 may contribute to tumorigenesis of HGSC through amplifying the master regulator MYC and increasing the dosages of these 46 genes related to OXPHOS.
DNA methylation is one of the most important epigenomic modifications, which undergoes extensive reprogramming during tumorigenesis and may impact gene expression (19, 20). We then analyzed the potentially regulatory effects of repeat element methylation, promoter methylation, and distal region methylation on gene expression.
We calculated the DNA methylation levels of whole genome and different annotated genomic elements (Fig. 3D; Supplementary Figs. S4B and S4C). For all patients, both satellites and long interspersed nuclear element 1 (LINE1) exhibited the most dramatic DNA demethylation (Fig. 3D), even for the cancer cells with global DNA methylation levels comparable with FTE (OC04, OC08, and lineage A1 in OC05; Supplementary Fig. S4B), indicating that strong demethylation of satellites and LINE1 is a hallmark of ovarian cancer. Next, to explore which DEGs were potentially regulated by DNA methylation of these genomic elements, we calculated the correlations between RNA expression and DNA methylation of these genomic elements (Supplementary Table S2). RNA expression of 61 DEGs upregulated in cancer cells were negatively correlated with DNA methylation levels of both satellites and LINE1 (Fig. 3E). Interestingly, the higher expression of genes associated with metallothioneins bind metals and interferon signaling, the two typical upregulated pathways in HGSC, was potentially affected by demethylation of satellites and LINE1 (Fig. 3E and G). In addition, we found that the upregulation of these 61 genes was associated with poorer progression-free survival (PFS) in HGSC (Fig. 3I). As for the DEGs downregulated in primary tumor cells, we revealed that the genes associated with translation elongation may result from demethylation of these two repetitive genomic elements (Fig. 3F and H); and the downregulation these genes predicted shortened PFS (Fig. 3I).
DNA methylation of gene promoters is one of the most important regulations of gene expression. Correlation and GO analysis show that genes associated with metallothioneins bind metals and interferon signaling are not only potentially regulated by satellite and LINE1 methylation, but also by their promoter methylation (Fig. 3J). Interestingly, the metabolic pathway of OXPHOS was significantly enriched but glycolysis was not (Fig. 3J), indicating that the expression of different metabolic pathways was regulated by different mechanisms.
The DNA methylation of distal regulatory regions, such as enhancers, may also be involved in gene expression regulation. Thus, we used MICMIC (21) to identify methylation regulation networks of distal regulatory regions. We totally detected 16,941 regulatory region-target pairs (Supplementary Table S3). Of these pairs, 53.9% were anticorrelated with the expression of target genes. Most target genes were potentially regulated by DNA methylation of less than five regulatory regions. These results demonstrate that some genes are not only potentially regulated by the methylation of their promoters, but also by methylation of distal regulatory regions. For example, the upregulated gene IFI27 (Fig. 3K), whose expression was negatively correlated with its promoter DNA methylation (Supplementary Table S3), was also negatively correlated with methylation of a distal region [∼235 kb away from IFI27 transcriptional start site (TSS); correlation, −0.26; P value, 4.8E−13]. In addition, this distal region is overlapped with a known ovarian enhancer (Fig. 3K) and its flanking regions exhibit higher chromatin accessibility in primary tumor cells compared with FTE cells, highlighting its putative regulatory roles of IFI27 expression.
Chromatin accessibility and the key transcription factors in HGSC
Binding of transcription factors (TF) to the cis-regulatory elements of a gene is a fundamental mechanism of transcription regulation, and TF access is often regulated by chromatin states. We found that the global chromatin accessibility of cancer cells tended to be more open than that of FTE cells (Fig. 4A), whereas the accessibility of regions adjacent to TSSs was less open in cancer cells than that of FTE cells (Fig. 4B). In addition, cancer cells had less and shorter nucleosome-depleted regions (NDR) than FTE cells (Fig. 4C), indicating that the cancer cells had less open chromatin states of focal open chromatin regions. The results of different genomic elements show that the accessibility of CpG islands (CGI), promoters, and microsatellite regions evidently decreased in cancer cells (Fig. 4D).
To explore potential key regulators in cancer cells, we performed TF binding motif enrichment of NDRs using chromVAR (22). In total, motifs of 107 TFs were more open in cancer cells, whereas 126 were more closed (Fig. 4E). To further explore which TFs were more likely to participate in DEG expression regulation, we used Lisa (23) to infer the transcriptional regulators of the DEGs (Fig. 4F). Of note, NFE2L2 was the most significantly enriched TF predicted using DEGs and simultaneously its binding motifs were more open in cancer cells than FTE evaluated by NDRs. Previous studies have verified that NFE2L2 can promote ovarian cancer growth (24) and it is associated with resistance to chemotherapy (25, 26). However, the expression level of NFE2L2 has no difference between cancer cells and FTE cells. Of these 107 TFs, expression levels of four TFs in cancer cells are significantly higher than that in FTE cells: FOXK1, TFAP2C, NR2F6, and DDIT3. FOXK1 is an important regulator to induce aerobic glycolysis (27, 28), which was previously proved to facilitate cell proliferation and metastasis in ovarian cancer (27, 29). Considering that the glycolysis pathway was upregulated in HGSC and the binding motifs of its transcriptional regulator FOXK1 exhibited high accessibility, we suggested that FOXK1 could be an attractive candidate therapeutic target for HGSC. For the 126 TFs with more closed chromatin states, several members of the GATA family were enriched. However, GATA1-5 were almost not expressed in cancer cells and FTE cells (Supplementary Fig. S5), whereas only GATA6 was expressed and its expression levels were different between cancer cells and FTE cells. Therefore, we considered that GATA6 played more important roles during tumorigenesis of HGSC among these GATA members. Previous studies had shown that loss of GATA6 expression led to nuclear deformation in ovarian cancer (30, 31). Here, we also illustrated that GATA6 may also regulate tumorigenesis process of HSGC through affecting the chromatin states of its target genes.
The gene expression profiles of abdominal metastases are similar to their matched primary tumors
After describing several key differences between FTE and primary tumors, we next explored the differences between primary tumors and matched metastases. On the basis of genetic lineage analyses, we compared the primary tumor cells and matched metastases within the same genetic lineage in the same patient.
Results show that very few genes were differentially expressed for the patients with abdominal metastases (Fig. 5A) and primary tumor cells exhibited close relationships to its matched metastatic cancer cells (Fig. 5B). These results propose that the microenvironment of abdominal metastases did not strongly affect the global RNA expression patterns of cancer cells, although the cancer cells have invaded and proliferated in a different organ, and imply that the metastatic capacity was probably acquired at an early stage during tumorigenesis. In contrast, previous studies based on bulk RNA profiling have reported several differences between primary tumors and metastases (32–34). Such contradiction may be caused by different technologies used in the researches. The DEGs between primary tumors and matched metastases revealed by bulk profiling may reflect the differences of different genetic lineages of cancer cells or different cell type compositions of the primary tumors and metastases, but not the actual differences of cancer cells from the same genetic lineage in the primary tumors and matched metastases.
Although only a few DEGs between metastases and primary tumors were shared by different patients (Fig. 5A), these genes may provide some clues to identify critical genes promoting or suppressing metastasis. For example, FN1 exhibited higher expression in metastasis of 2 patients (Fig. 5A). Previous studies have proved that FN1 can promote epithelial–mesenchymal transition and metastasis in serous ovarian cancer (35, 36). We next tested the functions of one of the five genes (HSPA6) downregulated in metastases shared by 2 patients (Fig. 5A). The results showed that overexpression of HSPA6 inhibited migration and invasion of the cancer cells, suggesting that its downregulation promotes tumor metastasis of HGSC (Fig. 5C–G).
Cancer cells maintained DNA methylation levels and patterns during abdominal metastasis
We found that global DNA methylation levels of cancer cells are stable during abdominal metastasis in most patients (Fig. 6A). Thus, we speculated that cancer cells maintained the global methylation characteristics of their genetic lineages during metastasis, so that the DNA methylation patterns can be used to describe cancer lineage histories. We then performed unsupervised hierarchical clustering of DNA methylation (see Supplementary Materials and Methods), and found that DNA methylation patterns produce similar lineage histories to those inferred by SCNAs (Fig. 6B; Supplementary Fig. S6). In contrast, different genetic lineages of cancer cells in the same patient could not be distinguished based on the clustering of the transcriptome and chromatin accessibility (Supplementary Fig. S6).
Intratumor heterogeneities of primary tumors provide insights into the metastasis process
Reconstruction of genetic lineages could unveil ITH and provide precise routes of metastasis. Different lineages not only have different SCNAs, DNA methylation levels, and RNA expression patterns but may also have different metastatic potential. For example, both patient OC09 and OC14 harbor two genetic lineages (lineage A1 and A2) of cancer cells, and the lineage A1 had established secondary tumors, but the other lineage (lineage A2) had not (Fig. 6A). Therefore, we wondered what kind of primary tumor cells possess greater metastatic potential. To explore the common upregulated genes in lineage A1 for patient OC09 and OC14, we performed differential expression analysis and revealed that CCN1 and HSP90AA1 were upregulated in the metastasized lineage A1 for both OC09 and OC14 (Fig. 7A), which prompted us to test if inhibition of CCN1 and HSP90AA1 can suppress the metastasis of ovarian cancer cells. We knocked out CCN1 using the CRISPR-Cas9 system in ovarian cancer cell line SKOV3 and ES2 (Fig. 7B; Supplementary Fig. S7B). Results showed that knockout of CCN1 significantly restrained the migration and invasion of cancer cells by wound healing assay and transwell assays (Fig. 7C and D; Supplementary Figs. S7A and S7C). We inhibited HSP90AA1 using the HSP90α inhibitor TAS-116. We revealed that TAS-116 can dose-dependently inhibit ovarian cancer cell migration and invasion (Fig. 7E–G; Supplementary Figs. S7D–S7F). In addition, TAS-116 could significantly suppress the growth of cancer cells in a time- and dose-dependent manner (Fig. 7H; Supplementary Fig. S7G). Together, our results demonstrate that CCN1 and HSP90AA1 may participate in the metastasis of HGSC and may serve as potential therapeutic targets for ovarian cancer.
To further unveil the shared upregulated pathways of the metastasized lineage between OC09 and OC14, we performed gene set enrichment analysis (GSEA; ref. 37). Using hallmark gene sets, we observed hallmarks of interferon gamma and alpha responses were significantly enriched in the metastasized lineage A1 (Fig. 8A–D). Consistently, results based on ontology gene sets show that metastasized lineage A1 exhibited higher expression associated with virus responses (Fig. 8E and F). In addition, lineage A1 also upregulated genes involved in blood vessel morphogenesis (Fig. 8E and F). In OC09, we further showed that the binding motifs of TFs associated with interferon signaling were more open in lineage A1 than A2 (Fig. 8G).
With respect to DNA methylation, we identified that, for both OC09 and OC14, the lineage with lower residual DNA methylation levels (lineage A2) did not metastasized. To further investigate which genomic regions were demethylated in nonmetastasized lineage A2, we performed enrichment analysis of the differentially methylated tiles between lineage A1 and A2. Results showed that the tiles with decreased DNA methylation in lineage A2 were strongly enriched in intergenic regions, satellites, and LTR (Fig. 8H).
Besides OC09 and OC14, primary tumors of four other patients possessed at least two genetic lineages, all of which metastasized or did not. Although these lineages present comparable metastatic potential, differential gene expression analysis shows that the phenotypes of these genetic lineages were distinct (Supplementary Figs. S8A–S8H).
Discussion
Here, we introduce an improved single-cell multiomics sequencing method, and provide the first systematic exploration about inter- and intratumor heterogeneities of SCNAs, DNA methylome, chromatin accessibility, and the transcriptome of primary tumors and matched metastases in HGSC at single-cell and single-base resolution.
Our single-cell multiomics sequencing technology exhibits several advantages in studying tumor biology. First, cell identities can be more accurately determined by multiomics data. Most studies distinguished malignant cells from nonmalignant cells using scRNA-seq data, mainly based on the gene expression of cell type-specific markers and RNA-inferred SCNAs (38, 39). In contrast, using the single-cell multiomics sequencing, we can not only determine the cell identities via RNA data but also profile higher resolution SCNA patterns by DNA sequencing data. Second, we can explore the dynamic changes of DNA methylome, chromatin states, and transcriptome during tumor development after reconstructing the genetic trajectories of cancer cells. Similarly, using the multiomics technology, Clark and colleagues utilized RNA data to reconstruct developmental trajectories of mouse embryonic stem cells and explored the epigenetic and transcriptomic properties of different cell states along these trajectories (40). Third, we can uncover potential inter-omic regulatory relationships, which are very difficult to identify by bulk sequencing or separate single-cell analyses of different molecular layers.
On the basis of genetic lineage analysis, we explored the dynamics of the epigenome and transcriptome during tumor promotion. We noticed that the metabolic pathways were highly reprogrammed in HGSC (Fig. 2B). A century ago, Warburg and colleagues proposed that tumors mainly metabolized glucose through glycolysis rather than OXPHOS (41). However, we showed that both OXPHOS and glycolysis were increased in ovarian cancer cells. Consistently, recent studies based on cell lines and mouse models have reported that ovarian cancer cells displayed higher rates of OXPHOS and glycolysis than normal epithelial cells (42, 43). We further revealed that these two pathways may be regulated by different mechanisms. The upregulation of genes related to OXPHOS may be due to either the Chr8 amplification to increase their dosages, or upregulate their upstream regulators (e.g., MYC), or their promoter demethylation; while the binding motif of the key TF involved in glycolysis, FOXK1, exhibited high accessibility in HGSC and may upregulate the glycolysis-related genes.
Taking advantage of single-cell multiomics sequencing technology, we revealed the high ITH in HGSC. We show that most tumors were composed of diverse lineages of cancer cells. These lineages vary in DNA methylation levels, RNA expression levels, and metastatic potential. The results of patient OC09 and OC14 show that the lineage with a higher residual DNA methylation level escaped from the primary tumor and invaded other tissues, whereas the other lineage, with a lower residual DNA methylation level, did not (Fig. 6A). Furthermore, the metastasized and probably more aggressive lineage exhibited upregulation of CCN1, HSP90AA1, and genes associated with virus response and angiogenesis (such as VEGFA). CCN1 has been proved to participate in cell proliferation, inflammation, and angiogenesis, and its overexpression was associated with poor prognosis of ovarian cancer (44, 45). Previous studies have shown that inhibition of HSP90AA1 may potentate the ability of chemotherapy (46–48). Here we revealed that CCN1 and HSP90AA1 may also be involved in ovarian cancer metastasis. While performing end point analysis to compare the primary tumors and metastases is a common means to investigate the mechanisms of metastasis, this approach cannot provide insights into the early events leading to metastasis. In our study, we used single-cell multiomics sequencing to compare primary tumor cells with different metastatic potential, which provided a new strategy for the study of the mechanisms of metastasis and help to better understand the metastasis process.
Next, we inspected the epigenetic and transcriptomic changes associated with tumor progression (Figs. 5 and 6). The results show that the global DNA methylation levels were maintained during metastasis. Unsupervised hierarchical clustering of DNA methylation shows that the DNA methylation patterns produce similar lineage histories to those deduced from SCNA patterns. As a covalent modification of genome, DNA methylation is relatively stable and involved in establishment and maintenance of cell type- and tissue-specific features during mitosis and development (49–51). Thus, despite the changes in DNA methylation at a relatively early stage of tumorigenesis, ovarian cancer cells retained the majority of their lineage-specific DNA methylation patterns during later tumor progression. The convergent pattern of DNA methylation and SCNAs was also demonstrated in prostate cancer (52). Another study has reported that phyloepigenetic patterns could recapitulate the phylogenetic histories reconstructed by SNVs (53). Therefore, the convergence of SCNAs, SNVs, and DNA methylome changes may exist in a variety of cancers.
In summary, using our single-cell multiomics sequencing, we obtained high-resolution profiles of HGSC, dissected its extensive inter- and intratumor heterogeneities, identified many potentially important variations during tumorigenesis, and explored the regulatory networks in HGSC. Our study provides highly valuable resources for understanding the molecular characteristics of HGSC, as well as for improving the diagnosis and individualized therapy of this disease.
Authors' Disclosures
K. Kee reports grants from Tsinghua University during the conduct of the study. No disclosures were reported by the other authors.
Authors' Contributions
Y. Wang: Conceptualization, resources, formal analysis, investigation, visualization, methodology, writing–original draft. H. Xie: Resources, data curation, software, formal analysis, visualization, methodology, writing–review and editing. X. Chang: Conceptualization, resources, formal analysis, supervision, funding acquisition, writing–review and editing. W. Hu: Formal analysis, investigation, writing–review and editing. M. Li: Formal analysis, investigation, methodology, writing–review and editing. Y. Li: Resources, investigation, methodology, writing–review and editing. H. Liu: Validation, visualization, methodology. H. Cheng: Investigation, methodology. S. Wang: Methodology. L. Zhou: Methodology. D. Shen: Methodology. S. Dou: Methodology. R. Ma: Methodology. Y. Mao: Methodology. H. Zhu: Methodology. X. Zhang: Methodology. Y. Zheng: Methodology. X. Ye: Methodology. L. Wen: Conceptualization, supervision. K. Kee: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, writing–review and editing. H. Cui: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, writing–review and editing. F. Tang: Conceptualization, resources, data curation, supervision, funding acquisition, validation, investigation, methodology, project administration, writing–review and editing.
Acknowledgments
The authors thank the Beijing Advanced Innovation Centre for Genomics, the Computing Platform of the Center for Life Science, and the National Center for Protein Sciences at Peking University for supporting our work to F. Tang. K. Kee and H. Cui were supported by the National Key Research and Development Program of China 2018YFA0107703 and 2016YFA0201400, respectively. X. Chang was supported by the National Natural Science Foundation of China (81971360 and 81671431).
The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.
Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).