RNA biomarkers from proximal liquid biopsy for diagnosis of ovarian cancer.

The RNA-sequencing dataset generated by the authors in the study is available at the European Genome-phenome Archive (EGA) under accession number EGAS00001005498.

Eighty-one utero-tubal lavage samples collected from 31 ovarian cancer patients and 50 patients with benign ovarian lesions were analyzed in this study. Upon RNA extraction and mRNA capture sequencing, 7 of the 81 samples (5 ovarian cancer and 2 control samples) were excluded for statistical analysis because of too low sequencing depth, resulting in a final cohort of 74 samples collected from 26 ovarian cancer patients with an average age of 61.5 years old and 48 patients with benign ovarian lesions with an average age of 63.8 years old. The demographic and clinical patient information is provided in Table 1 . Details for all samples with the reason of exclusion for further analysis is provided in supplemental Table 2. Table 1 Patient characteristics for utero-tubal lavage samples included in the transcriptomic analysis. Table 1 ovarian cancer ( n  = 26) control ( n  = 48) age in years, mean (min-max) 61.5 (46–78) 63.8 (51–83) stage  early stage (I-II) 4 -  late stage (III-IV) 22 - BRCA status  germline mutation 7 8  no mutation 17 0  unknown 2 40 indication for surgery  high grade ovarian cancer 26 -  benign ovarian mass/cyst - 9  menorrhagia - 3  pelvic organ prolapse - 6  leiomyomatous uterus - 2  normal endometrium - 6  RRBSO - 5  mature teratoma - 2  mucinous cystadenoma - 6  other - 9 Patient characteristics for utero-tubal lavage samples included in the transcriptomic analysis. Over all samples, the mapping rate to Ensembl genes was 90%, with a minimum of 20% and a maximum of 96% (Supplemental Fig. 1A). The total reads mapped to Ensembl genes varied from 1.2 million to 24.8 million per sample with a mean of 11.9 million reads per sample. Based on the coverage of artificial spike-in controls, the endogenous RNA mass per sample was calculated. The mean endogenous RNA mass detected per 1 mL fluid was 0.09 ng, with a minimum of 0.002 ng and a maximum of 0.72 ng. The endogenous RNA mass did not differ between the ovarian cancer group and the control group (Wilcoxon signed-rank test, two-sided, p  = 0.196, Supplemental Fig. 1B). Despite the high variability in mapping rate across the samples, RNA complexity of the samples was very stable. The total number of unique mRNAs ranged from 11,887 to 17,850 with a mean of 15,451 mRNAs per sample (Supplemental Fig. 1A). In total, 8139 genes were detected in all samples. To assess which tissues or cell types contribute mRNA molecules to the utero-tubal lavage fluid RNA profile, we evaluated tissue- and cell-type-specific mRNA signatures. The boxplots in Fig. 1 A highlight the relative contribution of tissues and cell types to utero-tubal lavage fluid compared to 23 human biofluids that were included in the Human Biofluid RNA Atlas [15] . Esophagus RNA markers were more abundant in utero-tubal lavage fluid than in the other biofluids, likely reflecting an epithelial RNA-signature shared between epithelial cells from esophagus and endometrium. Fallopian tube and ovary specific mRNAs were the second and third most enriched signatures, suggesting that the wash procedure enables detection of RNA originating from ovary and fallopian tube tissue. When comparing the relative RNA content of utero-tubal lavage to the 23 other biofluids, utero-tubal lavage fluid ranked as the twelfth highest concentrated fluid ( Fig. 1 B). The relative mRNA content of utero-tubal lavage was similar to that of ascites, broncho-alveolar lavage and platelet-rich-plasma. Utero-tubal lavage fluid contained 8-fold more RNA than platelet-free plasma, the most studied biofluid in the biomarker field. Fig. 1 Assessment of the tissues contributing RNA molecules to utero-tubal lavage fluid. (A) Boxplots showing the log2 fold change for a gene set with markers specific for a certain tissue or cell type. The log2 fold change is calculated between the median read count of all utero-tubal lavage samples and the median read count of all other biofluids. The tissues or cell types for which markers were selected based on the RNA Atlas Project are shown on the x -axis. (B) Barplot showing the relative mRNA content in 24 human biofluids. BAL, bronchoalveolar lavage fluid; CSF, cerebrospinal fluid; PFP, platelet-free plasma; PPP, platelet-poor plasma; PRP, platelet-rich plasma. Fig 1 Assessment of the tissues contributing RNA molecules to utero-tubal lavage fluid. (A) Boxplots showing the log2 fold change for a gene set with markers specific for a certain tissue or cell type. The log2 fold change is calculated between the median read count of all utero-tubal lavage samples and the median read count of all other biofluids. The tissues or cell types for which markers were selected based on the RNA Atlas Project are shown on the x -axis. (B) Barplot showing the relative mRNA content in 24 human biofluids. BAL, bronchoalveolar lavage fluid; CSF, cerebrospinal fluid; PFP, platelet-free plasma; PPP, platelet-poor plasma; PRP, platelet-rich plasma. Differential abundance analysis revealed 330 mRNAs that were significantly more abundant in utero-tubal lavage fluid from ovarian cancer patients compared to that from controls ( Fig. 2 A). Amongst the 330 mRNAs are bona-fide proliferation markers, such as Ki-67 and aurora kinase B (AURKB). A list with the results of the differential abundance analysis can be found in Supplemental Table 4. The normalized abundance of the 20 most differentially abundant mRNAs is shown in Fig. 2 B. Of note, principal component analysis of all expressed genes did not reveal clustering based on the clinical diagnosis of the donor (Supplemental Fig. 2). Gene set enrichment analysis of Hallmark gene sets demonstrated a statistically significant enrichment of four gene sets representing cell cycle deregulation: genes encoding cell cycle related targets of E2F transcription factors (normalized enrichment score (NES) = 2.34, padj = 0.02), genes involved in the G2/M checkpoint, as in progression through the cell division cycle (NES = 2.25, padj = 0.02), genes up-regulated by activation of hedgehog signaling (NES = 1.85, padj = 0.03) and genes that are important for mitotic spindle assembly (NES = 1.71, padj = 0.02; Fig. 3 ). PAX8, CA125 (MUC16) and HE4 (WFDC2), known lineage markers, were not differentially abundant between ovarian cancer patients and control samples (Supplemental Fig. 3). Differential abundance analysis on miRNA level revealed that 41 miRNAs were more abundant in the ovarian cancer group compared to the control group ( Fig. 2 B, Supplemental Table 4). Five of these miRNAs (let-7d-5p, miR-203a, miR-200b, miR-200c, miR-191) have previously been linked to the pathogenesis of ovarian cancer and were more abundant in plasma, serum or ascites of ovarian cancer patients compared to healthy controls [14] . Fig. 2 Differentially abundant mRNAs and miRNAs. (A) Volcano plot of differentially abundant mRNAs in ovarian cancer ( n  = 26) versus controls ( n  = 48). Upregulated genes with an adjusted p -value of less than 0.05 are shown in pink. No downregulated genes are detected. (B) Volcano plot of differentially abundant miRNAs in ovarian cancer ( n  = 26) versus controls ( n  = 48). Upregulated miRNAs with an adjusted p -value of less than 0.05 are shown in pink. No downregulated miRNAs are detected. (C) Boxplots comparing the Sequin spike normalized read counts per group for the top 20 most differentially abundant genes. The normalized read count per sample is shown as a dot. Samples obtained from ovarian cancer patients are pink, samples obtained from controls are blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article). Fig 2 Fig. 3 Enrichment plots of the four gene sets of the hallmark pathways that are enriched in ovarian cancer versus controls (adjusted p -value < 0.05). (A) Enrichment plot for genes encoding cell cycle related targets of E2F transcription factors (NES = 2.34, padj = 0.02). (B) Enrichment plot for genes involved in the G2/M checkpoint, as in progression through the cell division cycle (NES = 2.25, padj = 0.02). (C) Enrichment plot for genes upregulated by activation of hedgehog signaling (NES = 1.85, padj = 0.03). (D) Enrichment plot for genes important for mitotic spindle assembly (NES = 1.71, padj = 0.02). NES, normalized enrichment score; padj, adjusted p -value. Fig 3 Differentially abundant mRNAs and miRNAs. (A) Volcano plot of differentially abundant mRNAs in ovarian cancer ( n  = 26) versus controls ( n  = 48). Upregulated genes with an adjusted p -value of less than 0.05 are shown in pink. No downregulated genes are detected. (B) Volcano plot of differentially abundant miRNAs in ovarian cancer ( n  = 26) versus controls ( n  = 48). Upregulated miRNAs with an adjusted p -value of less than 0.05 are shown in pink. No downregulated miRNAs are detected. (C) Boxplots comparing the Sequin spike normalized read counts per group for the top 20 most differentially abundant genes. The normalized read count per sample is shown as a dot. Samples obtained from ovarian cancer patients are pink, samples obtained from controls are blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article). Enrichment plots of the four gene sets of the hallmark pathways that are enriched in ovarian cancer versus controls (adjusted p -value < 0.05). (A) Enrichment plot for genes encoding cell cycle related targets of E2F transcription factors (NES = 2.34, padj = 0.02). (B) Enrichment plot for genes involved in the G2/M checkpoint, as in progression through the cell division cycle (NES = 2.25, padj = 0.02). (C) Enrichment plot for genes upregulated by activation of hedgehog signaling (NES = 1.85, padj = 0.03). (D) Enrichment plot for genes important for mitotic spindle assembly (NES = 1.71, padj = 0.02). NES, normalized enrichment score; padj, adjusted p -value. Altered gene expression levels represent only a part of the complex transcriptional program in cancer cells. Alternative splicing, the differential inclusion and exclusion of exonic sequences in mRNA, is an additional mechanism that impacts the transcriptome. In a complementary analysis, we profiled differential exon usage in utero-tubal lavage fluid from cancer and control samples. Differential exon usage analysis revealed 407 exons that were significantly more abundant in utero-tubal lavage fluid from ovarian cancer patients compared to that from controls ( Fig. 4 ). A list with the results of the differential abundance analysis can be found in Supplemental Table 4. Of interest, 203 out of the 407 differential exons did not overlap with differentially abundant genes that were previously identified. Among these differentially abundant exons were exonic sequences belonging to TP53. TP53 is a tumor suppressor gene mutated in over 95% of all ovarian cancer cases, leading to either complete or partial loss of function. The exon segments that are differentially abundant in our cohort match with exon 5 of the main TP53 isoform [44] , which encodes for the highly conserved DNA-binding domain of the p53 protein and which contains the majority of the somatic mutations detected in ovarian cancer [45] . Also for MUC16, encoding the CA125 protein, and the oncogene Forkhead box M1 (FOXM1), increased abundance of selected exons was identified in samples from cancer patients. Fig. 4 Differentially abundant exons. (A) Volcano plot of differentially abundant exonic parts in ovarian cancer ( n  = 26) versus controls ( n  = 48). Upregulated exons with an adjusted p -value of less than 0.05 are shown in pink and yellow. Exonic sequences that belong to genes that are differentially abundant at gene level are shown in pink. Exonic sequences that do not belong to differentially abundant genes are shown in yellow. No downregulated genes are detected. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig 4 Differentially abundant exons. (A) Volcano plot of differentially abundant exonic parts in ovarian cancer ( n  = 26) versus controls ( n  = 48). Upregulated exons with an adjusted p -value of less than 0.05 are shown in pink and yellow. Exonic sequences that belong to genes that are differentially abundant at gene level are shown in pink. Exonic sequences that do not belong to differentially abundant genes are shown in yellow. No downregulated genes are detected. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fusion gene analysis of all 74 transcriptomes revealed a total of 414 raw fusion predictions in 26 ovarian cancer samples and 816 raw fusion predictions in 48 control samples. After stringent filtering, 64 high confidence fusion transcripts remained in the cancer group and 173 high confidence fusion transcripts were detected in the control group. A detailed list of the high confidence fusion transcripts is provided in Supplemental Table 5. The median number of high confidence fusion transcripts per control sample was 3 (min 0, max 11) and the median number of high confidence fusion transcripts per cancer sample was 2 (min 0, max 7). No significant difference in the number of fusion transcripts between both groups was detected (Mann-Whitney-U test, two-sided, p -value = 0.11). There was no overlap between the transcripts detected in our cohort and reported fusion transcripts in ovarian cancer [46] , [47] , [48] . Prior to the pre-processing step, the mRNA data contained 13,386 predictors, the miRNA data contained 1527 predictors and the exon data contained 509 predictors. The misclassification rate, sensitivity, and specificity for the individual data layers and for all possible combinations of the data layers are summarized in supplementary Table 6. The best classifiers were obtained with elastic net logistic regression when the 15 most significant markers were retained from each of the individual data layers after the pre-processing step. The multi-omics classifier based on combined mRNA, miRNA and exon data yielded the best performance, achieving an overall misclassification rate of 21%, a sensitivity of 66% and a specificity of 88%. A receiver operating characteristic (ROC) curve was built by repeating the analysis on 100 random splits of the data into training and test sets with a varying threshold ranging from 0 to 1 (Supplemental Fig. 4). The lowest misclassification rate was achieved with a threshold of 0.5 and resulted in an area-under-the-curve of 0.86. An overview of the most important predictors for the multi-omics classifier is provided in Supplemental Table 7.

E.H. is funded by the Fund for Scientific Research Flanders (FWO; 1133120N). Annelien Morlion is supported by the Special Research Fund (BOF) scholarship of Ghent University (BOF.DOC.2019.0047.01). This work is supported by Kom Op Tegen Kanker (Stand up to Cancer), the Flemish cancer society. K.L. is supported by the Israel Science Foundation (1104/17), the Israel Cancer Research Fund, Israel Cancer Association Early Detection grant, and the SPARK-Tel Aviv University grant.

Sample collection was approved by the ethics committee of Chaim Sheba Medical Center, Rabin Medical Center and Meir Medical Center, Israel (ClinicalTrials.gov identifier: NCT03150121 ). Written informed consent was obtained from each participant in accordance with the Helsinki declaration. Recruited patients underwent gynecological surgical procedures under general anesthesia, including hysteroscopy, hysterectomy and/or RRBSO. Eligible indications included high-grade ovarian cancer (primary or interval debulking), suspicious ovarian mass, risk reduction, or various other benign gynecological disorders. Utero-tubal lavage samples were collected before surgery, after induction of anesthesia, by surgeons in the three participating centers. An intrauterine insemination catheter (Insemi TM -Cath, Cook Inc. Bloomington, USA) or rigid pipelle uterine sampler (Endosampler, MedGyn, Addison, USA) was inserted into the endometrial cavity through the cervical canal. Ten ml of saline were flushed into the uterine cavity and fallopian tubes and immediately retrieved (at an average volume of 4.6 ml per patient). The utero-tubal lavage samples were immediately centrifuged at 480xg for 15 min to eliminate cells and the supernatants were stored at -80 °C. RNA extraction, library preparation and sequencing methods are described in detail in the supplemental methods section. Raw data is presented in supplemental Tables 1–7.

Utero-tubal liquid biopsy can be collected in a minimally invasive way and is an intriguing fluid to study in the context of ovarian cancer diagnosis, due to its contact with the epithelium of the fallopian tube, where these tumors arise. Here, we provide proof-of-principle that isolating RNA from utero-tubal lavage fluid is technically feasible. The mRNA capture sequencing data that was generated from utero-tubal lavage contains mRNA signatures specific for ovary and fallopian tube and can thus be used to explore liquid biopsy applications for ovarian cancer diagnosis. Bulk RNA sequencing allows to inspect RNA derived from the tumor as well as RNA representing the complex tumor-microenvironment. Our study revealed an upregulation of mRNAs involved in cell cycle regulation and proliferation in utero-tubal lavage fluid from ovarian cancer patients compared to control samples. Over 300 mRNAs were upregulated in ovarian cancer compared to control samples. V-Myb avian myeloblastosis viral oncogene homolog-like 2 (MYBL2), the most differentially abundant mRNA, showed an 8-fold upregulation in ovarian cancer patients compared with healthy donors. MYBL2 is a physiological regulator of cell proliferation, cell survival and cell differentiation, but it is frequently deregulated in ovarian cancer, contributing to tumorigenesis and progression [49 , 50] . The marker of proliferation Ki-67 is also among the upregulated genes, probably reflecting the persistent cell proliferation of ovarian cancer cells [51] . Two members of the E2F family of transcription factors (E2F1 and E2F8) are more abundant in ovarian cancer samples compared to controls. Deregulation of E2F transcription factors has been reported as a crucial player in ovarian cancer pathogenesis [51] , [52] , [53] and E2F1 has been suggested as therapeutic target. At exon level, centromere protein I (CENPI) is more abundant in the cancer group compared to the control group. CENPI is a known target gene of E2F1 that promotes chromosome instability in cancer [54] . Cell division cycle associated gene 5 (CDCA5) and aurora kinase B (AURKB) are upregulated in both utero-tubal lavage of ovarian cancer and in tumor tissue of ovarian cancer patients relative to benign ovarian tissue [55] . Both at gene and at exon level, no downregulated RNA markers were detected, which is in line with the hypothesis that ovarian cancer samples contain tumor derived RNA that is absent in the control samples. Beside identification of differentially abundant genes and exons, we also interrogated the presence of fusion genes in our dataset. In our cohort, high confidence fusion transcripts were detected in both cancer and control samples. Fusion gene analysis in ovarian cancer tissue and in ascites from relapsed patients has been reported in only a few studies and the contribution of fusions in this cancer entity remains unclear [46] , [47] , [48] . To our knowledge, this is the first time that RNA sequencing has been successfully applied to utero-tubal lavage samples of ovarian cancer patients to profile the extracellular RNA content. Barnabas et al. explored the proteomic profile of extracellular vesicles isolated from utero-tubal lavage fluid and constructed a 9-protein classifier for ovarian cancer diagnosis with 70% sensitivity and 76.2% specificity [25] . None of the 9 proteins that were included in the classifier, showed corresponding upregulation of mRNA in our cohort. It is known that the correlation between mRNA transcripts and generated protein expressions can be low due to differences in half lives and the post transcription machinery. Based on the available mRNA data, miRNA data and exon data, a multi-omics classifier was built to predict ovarian cancer. Combining the three different data layers resulted in the best classifier, with a sensitivity of 66% and specificity of 88%, indicating the added value of combining complementary data layers. A limitation of our study is that we did not have large enough cohort to segregate germline BRCA mutations carriers from BRCA -WT cases and controls, which could have highlighted a more robust classifier. This approach should be taken in future studies, since accumulating data hints that the expressional profile of BRCA -mutated müllerian epithelium significantly differs from the WT pattern. Another caveat to our study is that CA125 serum levels were lacking. As a result, the performance of our classifier could not be compared with the current gold-standard. A joint analysis of the transcriptomic and proteomic data could thus reveal useful insights that may not be deciphered from the separate analysis of mRNA or protein expressions. Biomarker development efforts to date clearly indicate that no individual biomarker, can provide sufficient sensitivity at high specificity for the early detection of ovarian cancer. In order to identify a robust multi-marker algorithm, it is necessary to explore alternative biofluids, such as utero-tubal lavage, and to combine different -omics approaches for biomarker discovery.

Ovarian cancer, the fifth leading cause of cancer-related mortality in women, a five-year survival rate below 45%, largely driven by late stage diagnoses [1] . In Europe, the incidence of ovarian cancer is 12.9 per 100,000 [2] . Ovarian cancer is often referred to as a ‘silent killer’ because local disease is usually asymptomatic and symptoms of advanced stage disease are nonspecific. More than 75% of affected women are diagnosed with metastatic disease that is rarely curable. Early detection of ovarian cancer is key as stage I disease has a 5-year survival rate of 93% [3] . Ovarian cancers are classified into histological subtypes, with various underlying transcriptional and mutational patterns. High-grade serous carcinoma is the most prevalent and most challenging subtype, which is common in genetically predisposed populations, such as germline BRCA1/2 mutation carriers, having an estimated lifetime risk of 54% and 23%, respectively. A risk-reducing bilateral salpingo-oophorectomy (RRBSO) around the age of 40, is the only effective approach to avoid ovarian cancer in these women, resulting in significant morbidity of early menopause [3 , 4] . Currently, no effective screening method for this cancer entity is available. Serum cancer antigen 125 (CA125), the most studied test for ovarian cancer screening, has important limitations [5 , 6] . Less than 50% of patients with early stage ovarian cancer have elevated CA125 levels, and elevated CA125 levels can also be observed in benign conditions, such as pelvic inflammatory disease, endometriosis, and ovarian cysts [7] . Human Epididymis Protein 4 (HE4) is another protein marker that has been evaluated as serum biomarker for ovarian cancer diagnosis and the test has received approval from the US Food and Drug Administration for use in women presenting with an ovarian mass in 2011 [8] . HE4 expression shows high specificity for ovarian cancer, however serum HE4 levels vary in smokers, in hormonal contraceptive users and the levels increase with aging [9] . Unfortunately, the use of CA125 and HE4 have not been effective in improving patient survival [10] , [11] , [12] , [13] . Therefore, there is a clinical need for non-invasive, robust, and reliable diagnostics for ovarian cancer detection. Extracellular RNAs (exRNAs) in blood and other biofluids have been identified as potential biomarkers for a wide range of diseases, including ovarian cancer [14] , [15] , [16] , [17] . These so-called ‘liquid biopsies’ may offer a non-invasive alternative to tissue biopsies for diagnosis, prognosis and treatment response monitoring. The biomarker potential of extracellular RNAs in serum, plasma, urine and ascites has previously been investigated for the early diagnosis of ovarian cancer. These studies mainly focused on selected microRNAs using reverse transcription quantitative polymerase chain reaction (RT-qPCR) [14] , while few applied RNA sequencing as a more unbiased and transcriptome wide approach [18] , [19] , [20] . There is increasing evidence that precursor lesions of high-grade serous carcinoma originate from the epithelium of the fallopian tube fimbriae rather than intraperitoneally. The fimbriae represent the distal end of the fallopian tube, adjacent to the ovaries [21] , [22] , [23] . The lag time from emergence of the first malignant cells to clinically overt high-grade ovarian cancer is approximately six years [21] , and shedding of tumor cells from ovarian cancer and its precursor lesions into the gynecological tract has been reported [24] . Sampling the cells of the fimbriae or their secreted biological products, through proximal liquid aspirated from the gynecological tract, may thus reveal markers of the initial lesions. Utero-tubal lavage fluid can be obtained after flushing saline into the uterine cavity and fallopian tubes and holds promise for a minimally invasive liquid biopsy technique as this can be performed during a routine office-visit at the gynecology department [25 , 26] . Previous studies looking into the biomarker potential of utero-tubal lavage fluid primarily focused on circulating mutant p53 DNA [26] and on proteomic profiling of extracellular vesicles isolated from utero-tubal lavage samples [25] . So far, the RNA content of this fluid remains to be investigated. The goal of this proof-of-concept study was to profile the extracellular transcriptome of utero-tubal lavage fluid using messenger RNA (mRNA) capture sequencing and small RNA sequencing to investigate the biomarker potential of extracellular mRNAs for ovarian cancer diagnosis.

The authors declare no potential conflicts of interest.