Transcriptome Profiling Reveals Matrisome Alteration as a Key Feature of Ovarian Cancer Progression

Background: Ovarian cancer is the most lethal gynecologic malignancy. There is a lack of comprehensive investigation of disease initiation and progression, including gene expression changes during early metastatic colonization. Methods: RNA-sequencing (RNA-seq) was done with matched primary tumors and fallopian tubes (n = 8 pairs) as well as matched metastatic and primary tumors (n = 11 pairs) from ovarian cancer patients. Since these are end point analyses, it was combined with RNA-seq using high-grade serous ovarian cancer cells seeded on an organotypic three-dimensional (3D) culture model of the omentum, mimicking early metastasis. This comprehensive approach revealed key changes in gene expression occurring in ovarian cancer initiation and metastasis, including early metastatic colonization. Results: 2987 genes were significantly deregulated in primary tumors compared to fallopian tubes, 845 genes were differentially expressed in metastasis compared to primary tumors and 304 genes were common to both. An assessment of patient metastasis and 3D omental culture model of early metastatic colonization revealed 144 common genes that were altered during early colonization and remain deregulated even in the fully developed metastasis. Deregulation of the matrisome was a key process in early and late metastasis. Conclusion: These findings will help in understanding the key pathways involved in ovarian cancer progression and eventually targeting those pathways for therapeutic interventions.


Introduction
Ovarian cancer (OC) is the deadliest gynecologic malignancy and a highly heterogenous disease [1][2][3]. High-grade serous OC (HGSOC) is the most prevalent and aggressive histologic subtype, constituting about 70% of all cases [4]. The other subtypes include endometroid, mucinous and clear cell OC. It was the first cancer to be extensively characterized by The Cancer Genome Atlas (TCGA) study [2], which found tumor protein p53 (TP53) mutations in almost all HGSOC tumors (96%). Among the other mutations, BRCA1 DNA repair associated (BRCA1) and BRCA2 DNA repair 2 of 22 associated (BRCA2) were found to be mutated in about 22% of samples, while seven other genes were mutated in 2-6% of tumors [2]. Another defining characteristic was the frequent copy number variations observed in HGSOC tumors. Based on the mRNA expression profiles, OC tumors were subclassified into differentiated, immunoreactive, mesenchymal or proliferative groups. Bowtell et al. had similarly identified six molecular subtypes of OC, which they termed C1-C6 [5]. They found that most high-grade tumors clustered into C1 (high stromal response), C2 (high immune signature), C4 (low stromal response) and C5 (mesenchymal, low immune signature). Together, these two studies provided a molecular overview of the aberrations at the genomic level as well as at the level of gene expression. Such molecular profiling provided a better perspective and opened up opportunities for personalized medicine. Integrated analyses of microarray based gene expression data sets have also resulted in the identification of prognostic signatures for OC [6]. However, these studies were limited to characterizing the primary tumors and did not directly reveal much about the progression and metastasis of OC.
For a more comprehensive understanding of the disease, one has to start from its tissue of origin. For many years, HGSOC was thought to originate from the ovarian surface epithelium. However, strong evidence for fallopian tube (FT) fimbria as the tissue of origin was uncovered when careful study of BRCA mutation carriers and later, HGSOC patients, revealed serous tubal intraepithelial carcinoma (STIC) in their FT fimbria as precursor lesions [7][8][9]. More precisely, a multi-platform genomic study revealed that even in those HGSOC cases, where STICs were not evident in the advanced stage, the potential site of origin remained the distal FT [10]. Even low-grade serous carcinoma of the ovary has now been reported to probably originate from the FT [11]. In recent years, significant progress has been made in our understanding of the genomic landscape of OC carcinogenesis through these sequencing studies. To a lesser extent, groups have also compared metastasis to primary tumors [12,13]. While the study by Grellety et al. was limited to comparing a panel of 429 genes in primary and metastatic tumors, Marchion et al. did a microarray analysis of ovarian surface epithelial cells, pelvic tumors and extrapelvic tumors. Another study based on RNA-seq of primary and metastatic tumors from omentum and bowel focused on phylogenetic analysis to determine the sequence of progression [14]. A study comparing matched primary tumors and omental metastasis for DNA copy number and mRNA expression changes identified the development of an aggressive phenotype with metastasis [15]. A recent report profiling gene expression before and after neoadjuvant chemotherapy showed marked changes in gene expression, including increased drug transport and peroxisomal pathways [16]. However, very few studies have done a comprehensive profiling of OC progression by comparing fallopian tube with primary tumor and metastasis, to identify the key pathways deregulated during carcinogenesis and subsequent progression. For example, a study based on exome sequencing of matched fallopian tube, ovarian tumor and metastasis from an OC patient revealed that transcoelomic metastasis did not result in much accumulation of genetic alterations [17]. In addition, a whole-exome sequencing of matched STICs, invasive fallopian tube carcinoma, ovarian tumors and omental metastasis identified diverse metastatic paths in OC, including STICs formed as a result of metastasis [18].
Due to these facts, a similar comprehensive approach comparing the transcriptome of fallopian tube with primary tumors and metastasis would be beneficial in revealing the oncogenic pathways altered, independent of genetic alterations. Since it has been demonstrated that accumulation of genetic changes is limited during metastasis [17], an RNA-seq based approach will be relevant to understand microenvironment induced gene expression changes. During metastasis, the OC cells encounter a new microenvironment that results in induction of adaptive changes [19]. Signals from the microenvironment have been shown to induce transcription factor and microRNA expression changes, which drive metastatic colonization, by changing the expression of their targets [20,21]. Similarly, interactions with the metastatic microenvironment can induce epigenetic changes such as DNA methylation through the induction of DNA methyl transferase [21]. While using RNA-seq for performing end point analysis of patient tumors is a relevant approach, it still does not provide any information of the intermediate steps and their regulation during metastasis. Since metastasis is a multi-step process, such information about the individual steps will be helpful in understanding and targeting the process that is mainly responsible for poor patient outcomes [22]. In vitro models, that can mimic certain steps of OC metastasis have been increasingly employed to understand the underlying mechanisms regulating these steps [23][24][25]. We have used such an organotypic three-dimensional (3D) culture model of the omentum-a common site of OC metastasis-to study the early regulation of metastatic colonization [20,21]. This model provides an insight into early metastatic colonization-the rate limiting step in metastasis [26]. In the present study, we have used a comprehensive approach of combining the endpoint RNA-seq analysis of OC patients' normal fallopian tubes, primary ovarian tumors and metastatic tumors with that of the transcriptomic changes occurring during early metastatic colonization, using the organotypic omentum culture model. This study not only reveals the key pathways deregulated during carcinogenesis and metastasis, but also provides a better understanding of the adaptive process induced during metastatic colonization. Such an approach revealed evidences that changes in the matrisome and related pathways are a defining feature of OC progression.

Study Design
RNA isolated from eight pairs of matched normal FT and primary tumors from OC patients were sequenced to identify the transcriptome changes during carcinogenesis. Similarly, RNA isolated from eleven pairs of matched primary and metastatic tumors were sequenced. Patient characteristics are included in Supplementary Table S1. This was combined with sequencing of RNA from three different HGSOC cells (Kuramochi, OVCAR4 and OVCAR8) seeded on the organotypic 3D culture model of the omentum [20,21] compared to controls, to identify the gene expression changes occurring during early metastatic colonization. Extensive analysis of the RNA-seq data, including differential expression, gene set enrichment, pathway analysis and interaction networks, was performed to identify the genes and pathways altered during carcinogenesis compared to metastasis as well as early metastatic colonization (3D omentum culture model) versus advanced metastasis (patient tumors).

Gene Expression Changes in Normal FT versus Primary Tumors and Primary Tumors versus Metastasis
A total of 2987 genes were differentially expressed in the primary tumor samples versus FT with a significance level of false discovery rate (FDR) < 5% and absolute fold change ≥ 2 ( Figure 1A Figure S1B). Butterfly plots demonstrating the positive and negative correlation of the genes in these two datasets are provided in Supplementary Figure S2. In order to identify the key deregulated genes, the top 50 upregulated and downregulated genes (fold change) were determined for primary tumor versus FT and for metastasis versus primary tumors. Heat maps were generated for these top 50 upregulated and downregulated genes ( Figure 1B).
Gene set enrichment analysis (GSEA) was done for the hallmark gene sets to identify the deregulated pathways during OC progression. 15 gene sets were positively correlated with primary tumors versus FT at FDR < 25% while four gene sets were negatively correlated at FDR < 25%. The most significantly correlated gene sets in primary tumors versus FT were E2f transcription factor (E2F) targets (enrichment score (ES) = 0.804), G2M Checkpoint (ES = 0.776) and MYC targets V1 (ES = 0.692) ( Figure 1C). Taken together, these results clearly indicate that during the process of carcinogenesis, the tumor acquires a proliferative phenotype where cell cycle regulation is affected along with the activation of oncogenic pathways.  Figure 1E), highlighting the importance of invasiveness and contractility during the process of metastasis. The E2F targets are important for carcinogenesis but are downregulated during metastasis (ES = −0.586) ( Figure 1F). Interestingly, G2M checkpoint (ES = −0.566) and MYC targets V1 (ES = −0.516) are inversely correlated to metastasis ( Figure 1F) showing that the primary tumors are more reliant on the proliferative phenotype.  Figure 1E), highlighting the importance of invasiveness and contractility during the process of metastasis. The E2F targets are important for carcinogenesis but are downregulated during metastasis (ES = −0.586) ( Figure 1F). Interestingly, G2M checkpoint (ES = −0.566) and MYC targets V1 (ES = −0.516) are inversely correlated to metastasis ( Figure 1F) showing that the primary tumors are more reliant on the proliferative phenotype. Other significant Hallmark gene sets positively correlating with primary tumors versus FT were mTOR complex 1 (MTORC1) signaling, MYC targets V2, DNA repair and oxidative phosphorylation (Supplementary Figure S3A). Additionally, apical junction, tumor necrosis factor a (TNFA) signaling via nuclear factor kappa B (NFKB), adipogenesis and angiogenesis (Supplementary Figure S3B) were negatively correlated. The metastasis versus primary tumors revealed additional gene sets that positively correlated with metastasis including apical junction, allograft rejection, interleukin 2 (IL2)-signal transducer and activator of transcription 5 (STAT5) signaling and adipogenesis (Supplementary Figure S3C). Similarly, those that were negatively correlated were MTORC1 signaling, protein secretion, androgen response and MYC targets (Supplementary Figure S3D).
The primary tumor versus FT and metastasis versus primary tumor data for the top deregulated genes ( Figure 1B) were further analyzed for pathways deregulated using Metascape [27] and presented as ontology clusters in Figure 2A,B, respectively. Muscle system process, core matrisome, blood vessel morphogenesis and extracellular matrix (ECM) organization were the most significant pathways in the primary tumors versus FT ( Figure 3A). Core matrisome, ECM organization, cilium movement and matrisome associated were the most significant pathways in metastasis versus primary tumors ( Figure 3B). Pathways involving matrisome, ECM organization, cellular signaling, development and morphogenesis were common in both of the data sets. In addition, an Ingenuity Pathway Analysis (IPA) was also conducted taking into consideration all the genes differentially expressed at 5% FDR with at least two-fold change in Primary tumors versus FT (  Figure S4A). The top canonical pathways included hepatic fibrosis/hepatic stellate cell activation, cAMP mediated signaling, G-protein coupled receptor signaling, calcium signaling and amyotrophic lateral sclerosis signaling. The top disease functions identified were cancer, cellular development, organismal injury and abnormalities (Table 1 (panel B)). The predicted upstream regulators included transforming growth factor beta 1 (TGFB1) (inhibited), vascular endothelial growth factor (VEGF), β-estradiol and erb-b2 receptor tyrosine kinase 2 (ERBB2) (activated) (  Figure  S4B). The top diseases and functions revealed by the analysis were digestive system development and function, embryonic development, and organismal development, which indicated that the process of metastasis mimicked regulatory mechanisms involved in development (Table 1 (panel E)). The key upstream regulators identified were myocardin (MYOCD) (activated), bone morphogenetic protein 4 (BMP4), RUNX family transcription factor 2 (RUNX2) and bone morphogenetic protein 2 (BMP2) (Table 1 (panel F)).

Overlapping Genes in Primary Tumors versus FT and Metastasis versus Primary Tumors
As the disease progresses from the formation of the primary tumor to the dissemination and development of the metastasis, it is expected to be associated with alterations in gene expression. Our analysis of sequencing data from patient specimens showed that the formation of primary tumors resulted in deregulation of 2987 mRNAs (Supplementary Table S2). The process of metastasis was accompanied by changes in 845 genes (Supplementary Table S3). We next proceeded to visualize the overlap in these two data sets. Figure 4A shows a volcano plot of the 2987 genes deregulated in the Primary tumors versus FT, where the genes that were also found to be differentially expressed in the metastasis versus primary tumors are highlighted as orange triangles. Similarly, the volcano plot of the 845 deregulated genes in metastasis versus primary tumors in Figure 4B highlights the genes that were also differentially expressed in Primary tumors versus FT as purple triangles. We then proceeded to visualize the extent of differential expression and the direction of change of the overlapping genes in both the datasets ( Figure 4C). The size of the dots in the volcano plot of the differentially expressed common genes in the metastasis versus primary tumors indicates the extent of differential expression of those genes in Primary tumors versus FT. Similarly, they are also color coded to indicate their upregulation or downregulation in FT versus primary ( Figure 4C). Of the 2987 and 845 differentially expressed mRNAs in Primary tumors versus FT and metastasis versus primary tumors, respectively, 304 mRNAs were overlapping ( Figure 4D, Supplementary Tables S4 and S5). Interestingly, an analysis of these overlapping genes revealed differences in the direction of change in the two data sets ( Figure 4E,F). While a majority of these genes were upregulated in metastasis versus primary tumors, most of them were downregulated in the primary tumors versus FT. This indicated that there were major differences in the processes involved in the initial carcinogenesis and the subsequent dissemination of the disease. The eight upregulated and 12 downregulated genes common to both the comparisons are listed in Table 2 (panels A and B). 231 genes were upregulated in metastasis versus primary tumors but downregulated in primary tumors versus FT. 53 genes were downregulated in metastasis versus primary tumors but upregulated in primary tumor versus FT. The top 10 oppositely regulated genes that are upregulated or downregulated in metastasis versus primary are listed in Table 2 (panels C and D), respectively.

Overlapping Genes in Primary Tumors versus FT and Metastasis versus Primary Tumors
As the disease progresses from the formation of the primary tumor to the dissemination and development of the metastasis, it is expected to be associated with alterations in gene expression. Our analysis of sequencing data from patient specimens showed that the formation of primary tumors resulted in deregulation of 2987 mRNAs (Supplementary Table S2). The process of metastasis was accompanied by changes in 845 genes (Supplementary Table S3). We next proceeded to visualize the overlap in these two data sets. Figure 4A shows a volcano plot of the 2987 genes deregulated in the Primary tumors versus FT, where the genes that were also found to be differentially expressed in the metastasis versus primary tumors are highlighted as orange triangles. Similarly, the volcano plot of the 845 deregulated genes in metastasis versus primary tumors in Figure 4B highlights the genes that were also differentially expressed in Primary tumors versus FT as purple triangles. We then proceeded to visualize the extent of differential expression and the direction of change of the overlapping genes in both the datasets ( Figure 4C). The size of the dots in the volcano plot of the differentially expressed common genes in the metastasis versus primary tumors indicates the extent of differential expression of those genes in Primary tumors versus FT. Similarly, they are also color coded to indicate their upregulation or downregulation in FT versus primary ( Figure 4C). Of the 2987 and 845 differentially expressed mRNAs in Primary tumors versus FT and metastasis versus primary tumors, respectively, 304 mRNAs were overlapping ( Figure 4D, Supplementary Tables S4 and S5). Interestingly, an analysis of these overlapping genes revealed differences in the direction of change in the two data sets ( Figure 4E,F). While a majority of these genes were upregulated in metastasis versus primary tumors, most of them were downregulated in the primary tumors versus FT. This indicated that there were major differences in the processes involved in the initial carcinogenesis and the subsequent dissemination of the disease. The eight upregulated and 12 downregulated genes common to both the comparisons are listed in Table 2 (panels A and B). 231 genes were upregulated in metastasis versus primary tumors but downregulated in primary tumors versus FT. 53 genes were downregulated in metastasis versus primary tumors but upregulated in primary tumor versus FT. The top 10 oppositely regulated genes that are upregulated or downregulated in metastasis versus primary are listed in Table 2 (panels C and D), respectively.
The common upregulated genes in metastasis and primary tumor include epiphycan (EPYC), collagen type X alpha chain 1 (COL10A1), fibronectin type III domain containing 1 (FNDC1) and podocan-like protein I (PODNL1), which are mainly extracellular matrix related genes. These findings further signify the importance of the tumor microenvironment in the progression of OC. The common downregulated genes in metastasis and primary tumor include anterior gradient protein 3 homolog (AGR3) and Sentan cilia apical structure protein (SNTN). Of the 2683 genes that were unique to primary versus FT and 541 genes unique to metastasis versus primary, the top 10 upregulated and downregulated ones are listed in Supplementary Table S6 (panels A-D). Thereafter, we proceeded to confirm the RNA-seq results by qRT-PCR (quantitative real-time polymerase chain reaction), for the top genes in each group of Table 2 in matched ovarian cancer patient primary tumors and metastasis ( Figure 5). The top three genes were chosen on the basis of the p-value and the published literature on their role in cancer. We also checked for their effect on ovarian cancer patient prognosis using KM Plotter (Figure 6). Patients were split by median expression of the gene and no exclusion criteria were applied (Supplementary Table S7). EPYC, COL10A1 and PLPP4, which were upregulated in both metastasis versus primary and primary versus FT were all confirmed by qRT-PCR. All three genes led to poor patient prognosis when overexpressed ( Figure 5A).    Among the genes downregulated in both the datasets, AGR3 and STN were confirmed to be downregulated in metastasis by qRT-PCR, but collagen type XXVIII alpha 1 chain (COL28A1) did not validate this ( Figure 5B). Similarly, decreased expression of AGR3 and STN resulted in poor patient prognosis, while COL28A1 had the opposite effect ( Figure 6B). Steroidogenic acute regulatory protein (STAR), nuclear receptor subfamily 0 group B member 1 (NR0B1) and serine protease 16 (PRSS16) are downregulated in metastasis versus primary but upregulated in primary versus FT. However, only STAR could be confirmed to be downregulated in metastasis by qRT-PCR ( Figure 5C). A decrease in STAR expression resulted in worse prognosis in patients, whereas NR0B1 and PRSS16 had moderate effects ( Figure 6C). Meis homeobox 3 (MEIS3), integrin subunit alpha 11 (ITGA11) and LYL1 basic helix-loop-helix family member (LYL1) are increased in metastasis and decreased in primary tumors versus FT. All of them were confirmed to be increased in metastasis by qRT-PCR ( Figure 5D). Among them, increased ITGA11 caused the worst prognosis ( Figure 6D). Among the genes downregulated in both the datasets, AGR3 and STN were confirmed to be downregulated in metastasis by qRT-PCR, but collagen type XXVIII alpha 1 chain (COL28A1) did not validate this ( Figure 5B). Similarly, decreased expression of AGR3 and STN resulted in poor patient prognosis, while COL28A1 had the opposite effect ( Figure 6B). Steroidogenic acute regulatory protein (STAR), nuclear receptor subfamily 0 group B member 1 (NR0B1) and serine protease 16 (PRSS16) are downregulated in metastasis versus primary but upregulated in primary versus FT. However, only STAR could be confirmed to be downregulated in metastasis by qRT-PCR ( Figure  5C). A decrease in STAR expression resulted in worse prognosis in patients, whereas NR0B1 and PRSS16 had moderate effects ( Figure 6C). Meis homeobox 3 (MEIS3), integrin subunit alpha 11 (ITGA11) and LYL1 basic helix-loop-helix family member (LYL1) are increased in metastasis and decreased in primary tumors versus FT. All of them were confirmed to be increased in metastasis by qRT-PCR ( Figure 5D). Among them, increased ITGA11 caused the worst prognosis ( Figure 6D).

Transcriptome Changes During Early Metastatic Colonization
The comparison of metastatic tumors with primary tumors of OC patients provides an end-point analysis of the differential gene expression occurring during metastasis. However, metastasis is a multi-step process and the knowledge of the regulation of the individual steps remains limited. The important, rate-limiting step of metastasis is the step of metastatic colonization [28]. At this step, the cancer cells reach the site of metastasis and then have to successfully adapt to the new microenvironment that they encounter [20][21][22]. Such adaptive processes involve productive interactions with the microenvironment and result in changes in gene expression, a process that is very poorly understood. Therefore, in vitro organotypic models, which mimic the microenvironment of the metastatic site, have been used to gain a better understanding of metastatic colonization [20,25,29,30]. Thus, we complimented our end point analysis of patient tumors with the gene expression analysis of HGSOC cells seeded on an organotypic 3D culture model of the human omentum ( Figure 7A) [20,21]. Green fluorescent protein (GFP) labeled HGSOC cells Kuramochi, OVACR4 and OVCAR8 were seeded on the organotypic 3D omentum culture model to mimic the early steps of metastatic colonization. This effectively replicates the initial interactions of the OC cells with the metastatic microenvironment and the resulting gene expression changes in the OC cells were analyzed by RNA-seq ( Figure 7A, Supplementary Figure S5). The 1182 differentially expressed genes represent the early changes induced by the new microenvironment of the metastatic site (Supplementary Table S8). The top deregulated genes, collagen type I alpha 2 chain (COL1A2), periostin (POSTN) and decorin (DCN) were then validated by qRT-PCR in all HGSOC cells ( Figure 7B). POSTN did not amplify in OVCAR4 control, but had an average Cq of 26.8 for OVCAR4 on 3D omentum culture, indicating a marked increase in expression. To identify the pathways involved in early metastatic colonization, a pathway analysis was performed using Metascape [27] for the differentially expressed genes in the HGSOC cells seeded on the organotypic 3D omentum culture model ( Figure 7C,E). The most significant pathways included ECM organization, integrin1 pathway, degradation of ECM, cell-cell adhesion via plasma membrane, collagen fibril organization and positive regulation of cellular component movement. Taken together, these indicate that the initial events during metastatic colonization involve interactions of the OC cells with the metastatic microenvironment, remodeling of the ECM and an invasive phenotype. Other pathways, such as the response to wounding, Ras, tumor necrosis factor (TNF) and platelet derived growth factor (PDGF) signaling, indicate early signs of reprogramming of the microenvironment and proliferative growth in close coordination with the microenvironment. The ontology clusters in Figure 7C help visualize the inter-relation between these pathways. While some of these pathways have a significant amount of cross-talk, others appear to constitute independent hubs. Some of these changes may be transient, while others may be sustained and essential for the metastatic tumor development. Therefore, this dataset was compared with the differentially expressed genes in the end-point analysis of metastasis versus primary tumors to identify those genes essential for both early colonization, as well as for the advanced metastasis ( Figure 7D). The 144 genes common to both the data sets represent those that are deregulated early during metastatic colonization, by the interactions with the metastatic microenvironment, and remain deregulated until the end point in the fully developed metastasis (Supplementary Table S9). These genes represent the essential drivers of the process of metastatic colonization and may serve as clinically relevant targets. The pathway analysis was also conducted for the 144 common differentially expressed genes in the metastasis versus primary tumors and HGSOC cells on 3D omentum culture ( Figure 7F). The most significant pathways identified were matrisome, core matrisome, ECM glycoproteins, extracellular matrix organization, matrisome associated, focal adhesion and integrin1 pathway. Apart from integrin signaling, PDGF signaling was also among the signaling pathways identified. Therefore, remodeling of the ECM, interactions with the microenvironment and alteration of the matrisome are the key features of the metastatic tumor development. These pathways also indicate the importance of the tumor microenvironment in the metastatic progression, and this needs to be considered for the development of new therapeutic strategies.

Discussion
In the current study, we did a comprehensive RNA-seq analysis of normal fallopian tubes, primary tumors and metastasis from OC patients, in combination with an organotypic 3D omental culture model for early metastatic colonization. While much information is available on the gene expression profiling of OC primary tumors compared to normal tissue, the information available about metastasis is still limited [31]. Therefore, a comprehensive study of carcinogenesis and progression together, to identify the common and unique signatures, is essential. Moreover, not much is known about the pathways deregulated during early metastatic colonization, the rate-limiting step of metastasis. Here, we combined end point analysis of primary versus metastatic tumors from patients with a 3D culture model of early metastatic colonization, to identify the common pathways deregulated during early colonization that remain relevant in advanced metastasis. Our results indicated that the core matrisome, extracellular matrix components, focal adhesion pathways, integrin mediated pathways and collagen formation pathways were clearly essential for this process.
It has been well documented that the genes coding for ECM components undergo deregulation in tumor progression [32][33][34][35]. Moreover, altered ECMs can drive tumor progression [36,37]. Tumorigenic properties of cancer cells such as migration and proliferation are largely influenced by the ECM composition and ultrastructure. Differential regulation of ECM components has been observed in the early and late metastasis stages, which further signifies the contribution of ECM in tumor progression [38]. Tumor extracellular matrix is contributed by both tumor cells and the stromal cells [38]. It has been demonstrated in OC that the dysregulation of ECM components and stiffness impacts the OC progression [39].
Collagen acts as a scaffold in tumor ECM and the synthesis and degradation of collagen largely impacts the metastatic attributes of cancer cells. Collagen is found in abundance in the ECM [40,41]. The deposition of collagen I, II and III is increased during cancer invasion and migration [42,43]. Rapid changes involving collagen deposition have also been observed in breast cancer. Collagen type I, II and IV have been utilized as potential biomarkers in colorectal cancer since degraded collagen fragments are released in the circulation [44]. Collagen VI has been found to be overexpressed in OC [45], while Collagen XI is greatly increased in the OC metastasis [46][47][48]. A comprehensive analysis of OC omental metastasis, including gene expression, matrisome, extracellular matrix organization, biomechanical properties, cytokine/chemokine levels and cellular profiles revealed matrisome changes are a defining feature and could be correlated with prognosis [49]. These findings also showed the close correlation between RNA-seq data with matrisome focused proteomics increase in ECM glycoproteins, such as fibronectin and fibrinogen, along with proteoglycans and secreted factors defined metastatic progression [49]. Our results indicate that core matrisome changes and ECM reorganization are key features of initial carcinogenesis and subsequent metastasis, including early metastasis (Figures 2  and 5). The remodeling of the tumor microenvironment was found to be an essential feature, consistent throughout the process of tumor initiation and progression.
Focal adhesions mediate the contact between cell and extracellular matrix, helping anchor the cell to the substratum and promote motility, playing a major role in metastasis [50][51][52][53]. Factors mediating or affecting the focal adhesions are usually upregulated in the tumors and this is also reflected in our sequencing data ( Figure 7F). The anchorage to ECMs is mediated by specific heterodimeric membrane glycoproteins-integrins, which interact with various ECM components such as fibronectin, collagen, vitronectin and laminin. Thus, elevated levels of specific integrins and their associated proteins are observed in the metastasis [54][55][56]. Tumor progression is also partly mediated by the vascular integrins, which are expressed in tumor vasculature [57,58].
Primary tumors differentially expressed 2987 genes as compared to the normal fallopian tubes. However, the process of metastasis resulted in changes in the expression of only 845 genes. Of these, 304 genes overlapped with those deregulated in primary tumors versus fallopian tubes. An interesting observation in our study was that these 304 common deregulated genes in primary tumors versus normal fallopian tube and metastasis versus primary tumors were not necessarily altered in the same direction ( Figure 4). This implied that certain changes necessary for the initial tumorigenesis were not essential for subsequent metastasis. Therefore, the process of cancer dissemination not only needed a new set of genes but the expression of certain genes that were overexpressed/repressed in the primary tumors had to be reversed. We went on to validate the top genes deregulated in similar and/or opposite directions by qRT-PCR in matched patient primary tumors and metastasis. Among the 12 genes tested, nine of them were validated by qRT-PCR ( Figure 5). Among them, increased expression of EPYC, COL10A1, PLPP4, MEIS3 and ITGA11 appeared clinically relevant as markers of poor prognosis (Hazard ratio 1.49, 1.49, 2.41, 1.47 and 1.68, respectively; Figure 6). On the other hand, increased expression of AGR3 was a positive prognosis marker ( Figure 6B). Of note, the datasets used for patient prognosis are mostly derived from primary tumor samples.
A key factor to consider is the different microenvironments of the primary tumor and the site of metastasis. Since gene expression changes can be regulated by the signals emanating from the microenvironment, this could be an important reason for the differences in metastasis. Such changes may not be reflected by genetic mutations in the cancer cells. The metastasizing cancer cells have to adapt to the new microenvironment and this is stimulated by the productive cross-talk with the microenvironment [28,59]. Our previous studies have shown that regulators of gene expression, such as microRNAs and transcription factors, are induced or repressed as a result of such reciprocal signaling [20,21]. Similarly, the cancer cells are also capable to induce mesothelial to mesenchymal transition in the mesothelial cells covering the omentum [24,60]. Our findings with the 3D omental culture model and comparison of patient primary and metastasis tumor samples, revealed that many of such initial adaptive changes continue to remain relevant in the advanced metastasis. These changes include multiple pathways involved in the interactions between the tumor cells with their microenvironment ( Figure 7F). These results support the recent observations made upon deconstructing omental metastasis [49], while demonstrating that many of these processes start really early during the process of metastatic colonization. These pathways are, therefore, needed for early colonization and remain essential for the subsequent progression into the advanced metastasis.

Reagents
Dulbecco's Modified Eagle Medium (DMEM), minimal essential medium vitamins, minimal essential medium nonessential amino acids, Trypsin and Penicillin-Streptomycin were obtained from Media Tech (Manassas, VA, USA).

Patient Samples
Matched pairs of fresh frozen primary tumors and normal fallopian tubes from eight OC patients were obtained from the Indiana University Simon Cancer Center Tissue Bank. Similarly, matched pairs of fresh frozen primary tumors and metastasis from 11 OC patients were also obtained from the Indiana University Simon Cancer Center Tissue Bank. RNA was isolated from these tissues using miRNeasy kit (Qiagen, Germantown, MD; Cat# 217004) following the manufacturer's protocol and submitted to the Center for Genomics and Bioinformatics core facility, Indiana University, Bloomington for library preparation and sequencing. The study was approved by the Institutional Regulatory Board of Indiana University (protocol numbers 1106005767 and 1606070934). Patient written consent was obtained and tissues were collected by Indiana University Simon Cancer Center Tissue Bank. Only deidentified patient specimens were provided by the tissue bank.

Cell Lines
Human HGSOC cell lines OVCAR8 were acquired from Ernst Lengyel at the University of Chicago and OVCAR4 was from Joanna Burdette, University of Illinois at Chicago. Kuramochi was procured from Japanese Collection of Research Bioresources. The cell lines used were genetically validated and tested to be mycoplasma free using respective services from Idexx BioResearch (Columbia, MO). The genetic validation was done using the CellCheck 16 (16 Marker STR Profile and Inter-species Contamination Test) and mycoplasma testing was done using Stat-Myco.

D Omental Culture and RNA Isolation
Total RNA was isolated from a panel of HGSOC cell lines (Kuramochi, OVCAR4 and OVCAR8) grown in 3D omental cultures as described previously [20]. Three independent experiments per HGSOC cell line were conducted. Briefly, human primary mesothelial cells (HPMCs) and normal omental fibroblasts (NOFs) were isolated from omentum of women donors undergoing surgery at Indiana University Health Bloomington Hospital. The NOFs and HPMCs were grown as primary cultures for 2-3 passages and used for the experiments. The NOFs (3.6 × 10 5 cells) were mixed with 91 µg Collagen I (BD Biosciences, San Jose, CA;Cat# 354236) and seeded in a 10 cm dish and allowed to attach for 4 h. Thereafter, 3.6 × 10 6 HPMCs were overlaid on them to form a confluent monolayer. This 3D omentum culture mimicking the surface layers of the omentum was allowed to grow for 24 h so that the cells could secrete ECMs and other factors to form a complex microenvironment. GFP labelled OC cells were then seeded on the 3D omentum culture to mimic the initial steps of colonization and allowed to grow for 48 h. The cancer cells were then isolated by fluorescence activated cell sorting and compared to control cancer cells directly seeded on 10 cm tissue culture dishes and sham sorted. Total RNA was isolated using miRNeasy kit (Qiagen, Germantown, MD; Cat# 217004) and submitted for RNA-seq. Only those genes that were significantly differentially expressed (< 5% FDR and absolute fold change ≥ 2) in all three HGSOC cells were considered for further analysis. We chose to use Kuramochi, OVCAR4 and OVCAR8 as they have been reported as 'likely/possibly' HGSOC cell lines [61]. Moreover, among the HGSOC cell lines, we have found them to be more suitable for functional assays and OVCAR8 also grows effectively in mouse xenografts [62,63]. Therefore, the same cells would be suitable for extending these studies in the future.

RNA Sequencing
Library preparation and next generation RNA sequencing was carried out at the Center for Genomics and Bioinformatics core facility, Indiana University, Bloomington. The library preparation was done using TruSeq Stranded mRNA HT Sample Prep kit (Illumina, San Diego, CA; cat#RS-122-2103) according to the manufacturer's protocol and 8-nucleotide barcodes were added for multiplexing. The barcoded libraries were cleaned by bead cut with AMPure XP beads (Beckman Coulter, Atlanta, GA; cat#A63882), verified using Qubit3 fluorometer (ThermoFisher Scientific, Waltham, MA) and 2200 TapeStation bioanalyzer (Agilent Technologies, Santa Clara, CA), and then pooled. The pool was sequenced on NextSeq 500 (Illumina) with NextSeq75 High Output v2 kit (Illumina, cat#FC-404-2005).

Analysis of Sequencing Data
Reads were adapter trimmed and quality filtered using Trimmomatic ver. 0.33 [64], with the cutoff threshold for average base quality score set at 20 over a sliding window of three bases. Reads shorter than 20 bases post-trimming were excluded. Cleaned reads mapped to human genome GRCh38 with gencode v25 annotation using Tophat2 ver 2.1.0 [65]. The number of reads mapped to annotated genes was counted using htseq-counts ver. 0.5.4p1 [66]. Significantly differentially expressed genes at 5% FDR with at least two-fold change were called using DESeq2 ver. 1.12.3 [67]. Categories based on overlaps among significantly differentially expressed genes between the two experiments (primary tumors versus FT and metastasis versus primary tumors) were analyzed using IPA (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis) and Metascape (metascape.org). Metascape is a web-based tool that has been under development since 2014 [27]. It enables multiple computational analysis of large scale data integrating gene annotation, membership analysis and has multi-gene-list meta-analysis capabilities. It allows the leveraging of 40 independent databases for studying functional enrichment, interactomes and gene annotation. Using Metascape, we ranked the deregulated pathways on the basis of -log 10 (p-values). Interaction modules were also clustered on the basis of functional similarity by Metascape.
Gene set enrichment analysis was conducted using GSEA Desktop ver. 3.0, build 0160 [68]. Hallmark gene sets from Molecular Signatures Database ver. MSigDB 6.1 were used. GSEArot algorithm [69] was used to compute significance of enrichment for each of the hallmark gene sets, taking into consideration that the primary tumors and the FT as well as the metastatic and primary tumors were matched specimens from individual patients.
Heat maps were generated using the heatmap.2 function from R programming tools (gplots). All RNA-seq data have been made available to the public (GEO accession numbers GSE137237, GSE137238, GSE137239).

Conclusions
Taken together, our RNA-seq data has generated a wealth of information in the form of crucial pathways and targets which can be further explored to unravel the molecular mechanisms regulating the OC metastasis. Overall, the results of this study emphasize the need to further explore the tumor microenvironment and extracellular matrix components in the progression of OC metastasis. OC patients typically present with the whole spectrum of metastasis ranging from early lesions to advanced metastasis. Targeting key components of the pathways identified in Figure 7F, or their regulators, could be an effective approach to treat the whole spectrum of metastasis.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-6694/11/10/1513/s1, Figure S1: Heat maps representing the deregulated genes in ovarian cancer, Figure S2: Butterfly plots, Figure S3: Gene set enrichment analysis (GSEA) for the Hallmark gene set in ovarian cancer, Figure S4: Ingenuity pathway analysis (IPA) of differentially regulated genes in ovarian cancer, Figure S5: Heat map representing differentially expressed genes in all 3 HGSOC cells (Kuramochi/OVCAR4/OVCAR8) seeded on the 3D omentum culture versus control. Table S1: Ovarian cancer patient characteristics and specimen information. Table S2: Differentially expresses genes in primary tumors versus FT. Table S3: Differentially expresses genes in metastasis versus primary tumors. Table S4: 304 differentially expressed genes in primary tumors versus FT that are also deregulated in metastasis versus primary tumors. Table S5: 304 differentially expressed genes in metastasis versus primary tumors that are also deregulated in primary tumors versus FT. Table S6: Common genes in metastasis vs. primary tumors and primary tumors vs. FT. Table S7: KM plotter settings for progression free survival analysis (PFS). Table S8: Differentially expressed genes in Kuramochi, OVCAR4 and OVCAR8 cells seeded on 3D omentum culture. Table S9: The 144 differentially expressed genes that are commonly deregulated in 3D omentum culture and in metastasis versus primary tumors.