Microarray expression profile and analysis of circular RNA regulatory network in malignant pleural effusion

ABSTRACT

Malignant pleural effusion (MPE) is a common complication of lung cancer. Accumulating evidence has suggested that circular RNAs (circRNAs) play important roles in oncogenesis and progression of cancer. However, the expression pattern of circRNAs in MPE remains largely unknown and awaits investigation. The study was designed to elucidate the potential roles of differentially expressed circRNAs in MPE. Herein, we detected a total of 1350 differentially expressed circRNAs and 1727 differentially expressed mRNAs in lung adenocarcinoma-associated malignant pleural effusion (LA-MPE) compared with tuberculous pleural effusion (TPE) by Clariom D Human Microarray. Among the top 5 up-regulated circRNAs (hsa_circ_0067705, hsa_circ_0025542, hsa_circ_0072793, hsa_circ_0084927, and hsa_circ_0085386), four were verified significantly up-regulated in LA-MPE by qRT-PCR and hsa_circ_0085386 had an increasing trend. CircRNA-miRNA-mRNA network for the top 5 up-regulated circRNAs was constructed and pathway analysis indicated that the enriched mRNA targets involved in PI3K-Akt signaling pathway, Axon guidance, Regulation of actin cytoskeleton and Rap1 signaling pathway were potentially regulated by these aberrantly expressed circRNAs. We generated specific circRNA profiles in LA-MPE for the first time. And analysis of circRNA regulatory network could provide evidence that circRNAs are important in MPE development because they participate in cancer-related pathways by sequestering miRNAs. Our findings suggested that aberrantly expressed circRNAs may be involved in the development of LA-MPE.

KEYWORDS:

• Malignant pleural effusion
• tuberculous pleural effusion
• circular RNAs
• circRNA-miRNA-mRNA network
• functional enrichment

Introduction

Malignant pleural effusion (MPE) is a common, debilitating complication of a number of cancers, and is defined as the accumulation of a great deal of exudates in the pleural cavity, accompanied by the presence of malignant cells. The majority of MPE is due to metastatic disease, most commonly lung cancer, especially lung adenocarcinoma [1]. The presence of MPE implies an advanced-stage disease with reduced life expectancy and quality [2,3], and it is associated with high morbidity and mortality [4]. Though the promising progress in cancer therapy, the treatment of MPE remains palliative, with median survival ranging from 3 to 12 months [5]. The high mortality of MPE patients is largely ascribed to a lack of reliable biomarkers for disease detection and treatment, and most importantly the underlying biology is poorly understood [6]. Pleural fluid cytologic examination is commonly used for the diagnosis of MPE. However, the mean sensitivity of this method remains low (only ~ 60% with a range of 40–87%), and its accuracy is affected by sample preparation, cytologist expertise, and tumor types [7]. Moreover, with China being heavily burdened with tuberculosis [8], it is sometimes difficult to differentiate MPE from tuberculous pleural effusion (TPE). The current methods for treating MPE, including needle aspiration, indwelling pleural catheters, and pleurodesis, are ineffective [1]. Hence, it is imperative to deepen the understanding of the pathogenesis of MPE and then identify effective diagnostic and therapeutic markers for MPE.

It is currently accepted that a combination of increased fluid production and impaired lymphatic outflow underlie the development of MPE, and considerable efforts have been invested in exploring the mechanism of tumor-host interplay [9] and T cell subsets [10–14] in MPE formation. The advent of next-generation sequencing technologies has fueled the discoveries of new classes of non-coding RNAs that may trigger the development of MPE research. For instance, microRNAs have been reported to be promising candidates for MPE biomarkers [15–18]. Recently, another class of non-coding RNA molecules, circRNAs, has drawn increasing interests for their regulatory roles in cancer and other diseases. CircRNAs are usually abundantly expressed, with tissue- and developmental stage-specific patterns [19–21] and the covalently closed continuous loop structure confers them the resistance to RNases [19,20,22], endowing them clear advantages as novel diagnostic markers in comparison to linear RNA [20,21]. In addition, it has been proposed that circRNAs can act as microRNA sponges, form RNA-protein complexes, and regulate target gene transcription and splicing [23,24]. All that matters is that circRNAs may play important roles in the carcinogenesis and progression of cancer, exhibiting distinctive expression profiles [25–28]. However, to date, the expression pattern and roles of circRNAs in MPE remain elusive.

In this study, we evaluated the expression profiles of circRNAs in MPE and TPE by Clariom D human microarray. We identified a number of differentially expressed circRNAs in MPE compared with TPE and constructed a predictive model for the circRNA-miRNA-mRNA regulatory network. We propose that aberrantly expressed circRNAs may be involved in the development of MPE by regulating target genes.

Materials and methods

Study subjects

The study protocol was approved by the Institutional Review Board for Human Studies, Beijing Luhe Hospital, Capital Medical University, Beijing, China. Informed consent was obtained from all subjects. A total of six pleural effusion (PE) patients were recruited for microarray assay. Three patients had MPE associated with recently diagnosed lung adenocarcinomata. MPE was diagnosed by the presence of malignant cells in pleural fluid and/or on a pleural biopsy specimen. The other three patients were diagnosed with TPE, as evidenced by the growth of Mycobacterium tuberculosis from PE or by demonstration of granulomatous pleuritis on a pleural biopsy specimen without evidence of other granulomatous diseases. After antituberculosis chemotherapy, disappearance of PE and clinical symptoms was observed in all patients with tuberculous pleurisy. At the time of sample collection, none of the patients had received any anticancer treatment, antituberculosis chemotherapy, corticosteroids, or other nonsteroid anti-inflammatory drugs. Patients were excluded if they had received any invasive procedures directed into the pleural cavity or if they had suffered chest trauma within 3 months prior to hospitalization. The demographic characteristics of the 6 patients are summarized in Table 1.

Table 1. Demographics of the study participants

Group	No.	Age (years)	Sex	BMI (Kg/cm^2)	Accompanying diseases	TP (g/L)	ALB (g/L)	ADA (U/L)	LDH (U/L)
	1	72	F	28	HP, CHD	49.6	35	16.1	315
LA-MPE	2	83	M	21.9	CI	48.0	27.5	25.7	1286
	3	81	M	28.4	HP	49.2	29	5.5	175
	1	53	F	25.1	HP	49.6	33.3	60.7	515
TPE	2	45	M	20.9	No	48.4	25.6	41.7	522
	3	85	F	16.9	CHD, AF	47.3	24.5	40.9	102

LA-MPE: lung adenocarcinoma associated malignant pleural effusion; TPE: tuberculous pleural effusion; HP: hypertension; CHD: coronary heart disease; CI: cerebral infarction; AF: atrial fibrillation; TP: Total protein; ALB: Albumin; ADA: adenosine deaminase; LDH: lactic dehydrogenase.

Sample collection and processing

The PE samples were collected from each subject in heparin-treated tubes, using a standard thoracocentesis technique within 24 hours after hospitalization. The obtained PE specimens were immediately immersed in ice and were then centrifuged at 1200 rpm for 5 mins. The cell-free supernatant of PE was immediately frozen at -80°C after centrifugation for subsequent determination of soluble mediators. The cell pellets were resuspended in TRIzol reagent (Invitrogen, NY, USA) to extract the total RNA according to the manufacturer’s instructions. RNA samples were stored at -80°C before use.

CircRNA, miRNA and mRNA microarray expression profiling

The purity and concentration of extracted RNA were determined by OD260/OD280 using a NanoDrop ND-1000 instrument (NanoDrop Thermo, Wilmington, DE, USA). The integrity of the RNA was assessed by 1% formaldehyde denaturing gel electrophoresis.

Microarray profiling was conducted by Premedical Technology Corporation, Beijing, China. Sample labeling, microarray hybridization and washing were performed according to the manufacturer’s standard protocols. Briefly, each purified RNA sample was transcribed to double-strand cDNA, followed by cRNA synthesis and biotin-labeling. The labeled cRNAs were hybridized onto the Clariom D Human Assay Microarrays, which included transcriptome-wide gene- and exon-level expression profiles. After washing, the arrays were scanned using the GeneChip Scanner 3000 7G (Affymetrix, Santa Clara, CA, USA). The Affymetrix GeneChip Operating Software was used to analyze the acquired array images. Quantile normalization and subsequent data processing were performed using the R software package (R version 3.4.4). Differentially expressed RNAs between MPE and TPE were estimated by fold-change filtering combined with Student’s t-test. False discovery rate (FDR, <0.05) was calculated in order to correct the P-value. Transcripts with fold change in expression ≥ 2.0 and a p-value < 0.05 were considered significantly differentially expressed. Hierarchical clustering was performed to show the distinguishable genes expression pattern among samples. Differentially expressed RNAs were filtered and illustrated as a volcano plot. Principal component analysis was performed for differentially expressed circRNAs and mRNAs obtained from the microarray assays.

Validation of candidate circRNAs using quantitative reverse transcription real-time polymerase chain reaction (qRT-PCR)

The validation of selected candidate circRNAs was performed using qRT-PCR. Total RNA was reversely transcribed to cDNA according to the manufacturer’s instructions. Quantitative RT-PCR analysis and data collection were performed on ABI 7900HT Real-time PCR System (Applied Biosystems). The primers were synthesized by Premedical Technology Corporation (Beijing, China) (Table 2). RPS18, GAPDH or 18SrRNA were used as the references for normalization. All PCRs were triplicated. The relative expressions of circRNAs were calculated using the formula 2^{−(ΔCt MPE – ΔCt TPE)}.

Table 2. Specific circRNA primers used for quantitative RT-PCR analysis

Name	Sequence
hsa_circ_0001554	F: AGTTGGGGTGGTTTGAGGTT
	R: GCCAAACTCAGGTTCATTTCCT
hsa_circ_0003954	F: TGAGGTCATGCAGTGTCACC
	R: TAGGGTCGCTCCATGGTTGG
hsa_circ_0004086	F: GTGACCCAGGTGGCACATTC
	R: AGTTGTTCATGCTCGGAGAGT
hsa_circ_0050579	F: TGGCCTCCAGGTTCCAATGA
	R: AAGCCTGCCAGGAACACAAG
hsa_circ_0067705	F: CCACTATGTCTTGATTCCCTCG R: CCACCAAGAGCCAAGAGGAT
hsa_circ_0025542	F: GCCCACCTGAATGACCTTG
	R: GGTAGATCCGTGCTCCGAC
hsa_circ_0085386	F: CGACATCCTGCGGACCTACTC
	R: GCTGCTTCCAATAAAAAGGCTC
hsa_circ_0072793	F: GAAGAAAGATGGACAGGTATGACA
	R: CCACCGCTGCTGTAACGA
hsa_circ_0084927	F: GCGTTGAGGAAGCATAAAGAC
	R: GTGAACCCCATGAGTACGAATA
18SrRNA	F: GTAACCCGTTGAACCCCATT
	R: CCATCCAATCGGTAGTAGCG
HS-RPS18	BN61101
HS-GAPDH	BN61104

Prediction of miRNAs and mRNAs for candidate circRNAs and computational analysis

Potential miRNA response elements (MREs) were searched for candidate circRNAs as miRNAs interact with circRNAs through MREs. MREs were predicted according to the cutoff value determined by the miRNA support vector regression (mirSVR) using miRanda [29]. The putative target miRNAs of the candidate circRNAs were identified by the aforementioned method. The target genes of miRNAs were further predicted based on TargetScan [30] and miRDB [31] to study the potential functions of the circRNAs. The predicted miRNAs and mRNAs were then compared with detected miRNAs and mRNAs as determined by the same microarray assay respectively. To study gene ontology (GO) enrichment, we used the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) to determine the top functional annotations and classifications of significantly differentially expressed circRNA-targeted gene pairs. In addition, Clusterprofiler package was used to annotate and categorize the functions of the target genes for the differentially expressed circRNAs in the pathways. The GO and pathway enrichment analyses were visualized by R software packages.

Construction of circRNA-miRNA-mRNA interaction networks

CircRNA was used as a bait to enrich the circRNA-miRNA-mRNA network according to the target prediction software miRanda [29], TargetScan [30] and miRDB [31]. To visualize their interactions, Cytoscape software was used to build the circRNA-miRNA-mRNA networks.

Statistical analysis

Unpaired Student’s t-tests were evaluated for significantly differentially expressed transcripts between MPE and TPE. Evaluation of statistical significance for the qRT-PCR results was performed by Student’s t-test using the SPSS Statistics 19.0 software package (IBM, Chicago, IL, USA). p < 0.05 was considered statistically significant. Plotting of the data was performed with Cytoscape, GraphPad Prism 5.0 (GraphPad Software, La Jolla, CA), or R software package (R version 3.4.4).

Results

Demographic characteristics of the study subjects

We recruited 6 pleural effusion patients to perform microarray assay, including 3 lung adenocarcinoma associated malignant pleural effusion (LA-MPE) and 3 tuberculous pleural effusion (TPE) at Beijing Luhe Hospital, Capital Medical University, Beijing, China. The demographic characteristics of the 6 patients are summarized in Table 1. All the patients were diagnosed based on histopathological examination. Most (5/6) of the patients had a history of accompanying diseases that had been stabilized. Both the LA-MPE and TPE are exudates on the basis of Light’s criteria [32]. The concentration of adenosine deaminase (ADA) in LA-MPE is much lower than that in TPE (p < 0.05).

Profiles of differentially expressed circRNAs and mRNAs in MPE and coding genes-mediated GO and pathway analysis

The Clariom D Human Array was used to profile circRNA and mRNA expressions for LA-MPE and TPE samples. In total, we detected 12,960 circRNAs (8,587 up-regulated, 4,373 down-regulated) and 20,666 mRNAs (13,343 up-regulated, 7,323 down-regulated). Among these, 1,350 circRNAs (817 up-regulated, 533 down-regulated) and 1,727 mRNAs (1,104 up-regulated, 623 down-regulated) were significantly differentially expressed (fold change ≥ 2.0 and p value < 0.05). Hierarchical clustering analysis on these circRNAs and mRNAs showed distinctive expression patterns for LA-MPE and TPE samples (Figure 1(a,b)), suggesting the regulation of LA-MPE pathogenesis is distinct from TPE. The differentially expressed circRNAs and mRNAs between the two groups were identified through using volcano plot filtering (Figure 1(c,d)). In both differentially expressed circRNAs and mRNAs, a majority of the transcripts (> 60%) were up-regulated in LA-MPE (Figure 2(a)). The differentially expressed circRNAs are widely distributed on chromosomes, including sex chromosome (Figure 2(b)).

Figure 1. Differential expression of circRNAs and mRNAs in MPE and TPE. (a, b) Hierarchical clustering analysis of differentially expressed circRNAs (a) and mRNAs (b) between MPE group and TPE group; each group included three individuals (FC (absolute value) ≥ 2.0, P < 0.05). Each column represents one sample, and each row represents a transcript. “Red” indicates high relative expression, and “Green” indicates low relative expression. (c, d) Volcano plots showing differentially expressed circRNAs (c) and mRNAs (d) in MPE. Red and green points represent up-regulated and down-regulated transcripts in MPE, respectively (FC (absolute value) ≥ 2.0, P < 0.05). MPE, Malignant Pleural Effusion; TPE, Tuberculous Pleural Effusion.

Figure 2. Characterizations of circRNAs and mRNAs expression in MPE compared with TPE. (a) Pie chart of aberrantly expressed circRNAs and mRNAs in MPE. (b) The distribution of differentially expressed circRNAs in human chromosomes. MPE, Malignant Pleural Effusion; TPE, Tuberculous Pleural Effusion.

We performed gene ontology analysis for differentially expressed mRNAs. As for up-regulated mRNAs in MPE (Figure 3(a)), the biological process analysis indicated that they were mainly involved in cell migration, extracellular matrix organization, cell-cell adhesion, and cell adhesion. The cellular component analysis showed that they were mainly involved in extracellular exosome, focal adhesion, cell surface, plasma membrane, and extracellular space. While the molecular function revealed that they mainly took part in cadherin binding involved in cell-cell adhesion, integrin binding, and protein binding. On the other hand, the GO enrichment analysis on the down-regulated mRNAs of MPE (Figure 3(b)) showed that the biological process was primarily enriched in T cell activation and regulation of immune response. The major cellular components were plasma membrane, T cell receptor complex, immunological synapse, and cytoplasm. The molecular function enriched in SH3/SH2 adaptor activity.

Figure 3. GO and KEGG analysis of differentially expressed mRNAs in MPE. (a, b) GO annotation of up-regulated (a) and down-regulated (b) mRNAs. BP, biological process; CC, cellular component; MF, molecular function. (c, d) KEGG pathway enrichment analysis for up-regulated (c) and down-regulated (d) mRNAs. The related pathways were ranked according to GeneRatio (Selection Counts/Selection Size).

KEGG analysis showed that there were 26 pathways related to up-regulated mRNAs in MPE, including PI3K-Akt signaling pathway, proteoglycans in cancer, focal adhesion, and Rap1 signaling pathway. The related pathways were ranked according to GeneRatio. The top 3 related pathways were PI3K-Akt signaling pathway, proteoglycans in cancer and focal adhesion (Figure 3(c)). Similarly, the down-regulated mRNAs were enriched in 37 pathways, such as Th17 cell differentiation, T cell receptor signaling pathway, Th1 and Th2 cell differentiation, Natural killer cell mediated cytotoxicity and so on. The top 3 related pathways were Cytokine-cytokine receptor interaction, Th17 cell differentiation and Natural killer cell mediated cytotoxicity according to GeneRatio (Figure 3(d)).

Robustness of microarray assay and validation of the differentially expressed circRNA candidates

To confirm the validity of microarray assay, we evaluated whether the identified differentially expressed circRNAs and mRNAs can distinguish MPE from TPE. We performed principal component analysis (PCA) on the differential circRNAs obtained from microarray data and confirmed that MPE and TPE were divided into distinctive groups (Figure 4(a)). Similarly, the PCA plot of the differentially expressed mRNAs showed that MPE were also separated from TPE (Figure 4(b)). Hence, the differentially expressed circRNAs and mRNAs could distinguish MPE from TPE, firmly supporting the robustness of microarray assay for discovering specific circRNAs for MPE.

Figure 4. The principle component analysis (PCA) on the microarray data. (a) PCA was performed on the differentially expressed circRNAs obtained from microarray data. The MPE samples can be separated from the TPE samples. (b) PCA plot of the differentially expressed mRNAs.

To confirm the microarray results, we firstly randomly selected four circRNAs (hsa_circ_0003954, hsa_circ_0004086, hsa_circ_0050579, and hsa_circ_0001554, FC ≥ 2.0, P < 0.05) for qRT-PCR validation. The primers used for the qRT-PCR analysis are listed in Table 2. All samples were evaluated in triplicates with negative and positive controls. Consistent with the circRNA microarray result, two circRNAs (hsa-circ-0001554 and hsa-circ-0004086) were confirmed to be up-regulated in LA-MPE (Figure 5(a)). Hsa-circ-0003954 and hsa-circ-0050579 failed qRT-PCR analysis with non-specific amplification. To further confirm the reliability of the microarray expression data and exclude technical reasons, we used an independent set of samples consisting of 9 LA-MPE and 9 TPE to verify the top 5 (relatively high fold change from the microarray assay) significantly up-regulated circRNAs (hsa_circ_0067705, hsa_circ_0025542, hsa_circ_0072793, hsa_circ_0084927 and hsa_circ_0085386). As shown in Figures 5(b) and (c), the relative expression of hsa_circ_0067705 and hsa_circ_0025542 was much higher in MPE than that in TPE, respectively. The relative expression of hsa_circ_0085386 in MPE had an increasing trend compared with that in TPE, though not statistically significant (Figure 5(d)). Of note, hsa_circ_0072793 and hsa_circ_0084927 were detected in MPE group, but not in TPE group (Figure 5(e,f)). Therefore, the top 5 significantly up-regulated circRNAs were selected for further analysis.

Figure 5. Validation of differentially expressed candidate circRNAs. The expression levels of circRNAs were detected using qRT-PCR. (a) Validation of the expression of has_circ_0001554 and has_circ_0004086. The heights of the columns represent the fold change (MPE/TPE) computed from the microarray data and the qPCR. (b, c, d, e, f) Validation of the expression of hsa_circ_0067705, hsa_circ_0025542, hsa_circ_0085386, hsa_circ_0072793, and hsa_circ_0084927. CircRNA expression was quantified using the 2^{−(ΔCt MPE – ΔCt TPE)} method with normalization to 18SrRNA expression levels. Values are expressed as the mean ± SEM (n = 9). qRT-PCR, quantitative reverse transcription real-time polymerase chain reaction; SEM, standard error of the mean; MPE, Malignant Pleural Effusion; TPE, Tuberculous Pleural Effusion; ND, non-detected.

Prediction of the top 5 candidate circRNA targeted genes and bioinformatics analysis for the interaction network of circRNA-miRNA-mRNA

To predict the targets of the top 5 significantly up-regulated circRNAs (hsa_circ_0067705, hsa_circ_0025542, hsa_circ_0072793, hsa_circ_0084927 and hsa_circ_0085386), we firstly established the basic circRNA-miRNA connectivity with miRanda. As mentioned above, the majority of circRNAs and mRNAs were up-regulated, and pathway analysis showed up-regulated mRNAs were mainly enriched in cancer-related pathways. To construct competing endogenous RNA (ceRNA) network, the predicted miRNAs were subjected to the down-regulated miRNA expression filter using data from the same microarray. After being filtered, the intersection miRNAs were ranked according to the mirSVR scores, thereby identifying the 5 highest ranking miRNAs for further analysis. There were only three miRNA candidates for hsa_circ_0025542, namely, hsa-miR-298, hsa-miR-607, hsa-miR-4710. The 5 highest-ranking miRNA candidates binding targets of the other four circRNAs were listed below: 1) For hsa_circ_0067705: hsa-miR-298, hsa-miR-3185, hsa-miR-544a, hsa-miR-4675, and hsa-miR-6071; 2) For hsa_circ_0072793: hsa-miR-4432, hsa-miR-4522, hsa-miR-1261, hsa-miR-6126, and hsa-miR-3123; 3) For hsa_circ_0084927: hsa-miR-4319, hsa-miR-520e, hsa-miR-520b, hsa-miR-623, and hsa-miR-302e; 4) For hsa_circ_0085386: hsa-miR-4692, hsa-miR-3612, hsa-miR-4698, hsa-miR-4514, and hsa-miR-5695. Following miRNA target predictions, we proceeded to generate the miRNA-mRNA interaction network with TargetScan and miRDB. After cross validation with the up-regulated mRNA expression data from the same microarray, we finally obtained a total of 23 down-regulated miRNAs and 675 up-regulated mRNAs that were predicted to be regulated by the top 5 up-regulated candidate circRNAs. To gain insights into the functions of the top 5 up-regulated candidate circRNAs in LA-MPE, we conducted Gene Ontology and pathway enrichment analyses for the predicted target genes. The Gene Ontology analysis showed that the target genes were mainly enriched for the GO terms of protein binding, plasma membrane, and extracellular exosome (Figure 6(a)). The pathway analysis revealed that the target genes were mainly enriched in PI3K-Akt signaling pathway, followed by Axon guidance, Regulation of actin cytoskeleton, and Rap1 signaling pathway (Figure 6(b)). Finally, based on the predictions of the miRNA and mRNA targets for the top 5 up-regulated candidate circRNAs and microarray data, we constructed a ceRNA network of circRNA-miRNA-mRNA interactions with Cytoscape software delineating the putative regulations on these targets in the pathogenesis of LA-MPE (Table 3 and Figure 7).

Table 3. The top 5 up-regulated circRNAs targeted circRNA-miRNA-mRNA network

CircRNAs	MRE1	MRE2	MRE3	MRE4	MRE5	No. of mRNA targets
hsa_circ_0067705	hsa-miR-298	hsa-miR-3185	hsa-miR-544a	hsa-miR-4675	hsa-miR-6071	142
hsa_circ_0025542	hsa-miR-298	hsa-miR-607	hsa-miR-4710			163
hsa_circ_0072793	hsa-miR-4432	hsa-miR-4522	hsa-miR-1261	hsa-miR-6126	hsa-miR-3123	154
hsa_circ_0084927	hsa-miR-4319	hsa-miR-520e	hsa-miR-520b	hsa-miR-623	hsa-miR-302e	36
hsa_circ_0085386	hsa-miR-4692	hsa-miR-3612	hsa-miR-4698	hsa-miR-4514	hsa-miR-5695	180

MRE: miRNA response element.

Figure 6. GO and KEGG analysis of the top 5 up-regulated circRNA targeted mRNAs. (a) GO enrichment for targeted mRNAs. All significant enriched GO terms were listed. (b) KEGG pathway enrichment for targeted mRNAs. Top 10 significant enriched pathway terms were listed.

Figure 7. The circRNA-miRNA-mRNA network. 23 miRNA and 675 genes were targeted by the top 5 significantly up-regulated circRNAs. Diamond, circRNA; Rectangle, miRNA; Ellipse, mRNA; Red, up-regulated, Fold Change (FC) ≥ 2.0 and p-value < 0.05; Green, down-regulated, FC≤-2.0 and p-value < 0.05.

Discussion

Malignant pleural effusion is an advanced-stage manifestation of lung cancer, especially lung adenocarcinoma [1,2]. Despite its longstanding recognition, there is currently no effective treatment option beyond needle aspiration to relieve the symptom of shortness of breath. Displaying versatility in regulating many pathological processes, non-coding RNAs may serve as another layer of mechanisms for the pathogenesis of MPE and be promising targets for therapy. In the current study, we profiled differentially expressed circRNAs in LA-MPE compared with TPE and constructed a regulatory network for circRNA-miRNA-mRNA in LA-MPE. We have identified specific circRNAs that may regulate the LA-MPE pathogenesis.

Non-coding RNAs play multidimensional roles in regulating many physiopathological processes [33]. Recently, many studies explored aberrantly expressed miRNAs in MPE and proposed their potentials as markers for diagnosis and/or prognostic predictions [15–18,34]. CircRNAs possess covalently closed loop structures and lack free 3ʹ- and 5ʹ- ends, therefore are more stable than other RNA molecules. This makes them more attractive candidates for biomarkers than other types of RNA. To date, the expression pattern and role of circRNAs in the MPE pathogenesis remains largely unknown. Some studies investigated the expression profiles of circRNAs in peripheral blood to extrapolate their potential roles in certain diseases [35–38]. But the conclusions should be carefully extended to tissues where the diseases reside, because expression of circRNAs is usually tissue-specific. To explore the role of circRNAs in MPE, we profiled circRNAs expression in LA-MPE and TPE clinical samples with microarrays. We identified 1,350 significantly differentially expressed circRNAs in LA-MPE, suggesting circRNAs may be important in regulating the pathogenesis of LA-MPE. What was noteworthy was that a majority of differentially expressed circRNAs were up-regulated. This result was in accordance with the clustering results on the differentially expressed mRNAs, where most of the genes were up-regulated in LA-MPE (Figure 2(a)). The fact that most genes were up-regulated in LA-MPE suggests that numerous pathways may promote the development and progression of MPE. Indeed, KEGG pathway analysis on differentially expressed up-regulated mRNAs confirmed that cancer-related pathways of PI3K-Akt signaling pathway, Proteoglycans in cancer, Focal adhesion, and Rap1 signaling pathway were enriched in LA-MPE.

The aberrantly expressed circRNAs may regulate the expression of parental genes and play important roles in the pathological processes of diseases. And a considerable number of studies have indicated that circRNAs can act as “miRNA sponges” to regulate gene expression [39–41]. As reported by the present study, most circRNAs and mRNAs were up-regulated, we speculated that circRNAs may regulate the expression of target genes by ceRNA mechanism. Based on the miRNA target predication software, we initially predicted miRNA targets of the top 5 up-regulated and validated circRNAs and then subjected them to down-regulated MPE-associated miRNAs，indicating up-regulated circRNAs might inhibit the expression of down-regulated miRNAs. Among the 23 potential miRNAs, some have been reported as important regulaters in many kinds of cancers [42–48]. For example, L Bao et al revealed that down-regulation of miR-298 increased P-Glycoprotein expression and induced doxorubicin resistance in sensitive breast cancer cells, suggesting miR-298 was associated with the chemoresistant mechanisms of metastatic human breast cancer by directly modulating P-Glycoprotein expression [42]. Similarly, a study by S Wei et al showed that miR-623 was down-regulated in lung adenocarcinoma and suppressed the invasion and metastasis targeting Ku80 through ERK/JNK inactivation mediated downregulation of MMP-2/9 [48]. No functional data is yet available for these miRNAs in MPE and may be considered promising candidates for future study. Subsequently, to further evaluate the biological roles of the top 5 up-regulated circRNAs, we predicted the target genes of these 23 down-regulated miRNAs. After subjecting predicted mRNAs to up-regulated LA-MPE associated mRNAs from the same microarray, we finally constructed a circRNA-miRNA-mRNA network consisting of 5 up-regulated circRNAs, 23 down-regulated miRNAs and 675 up-regulated mRNAs (Figure 7).

Functional annotation of target genes can provide insights of the biological functions of the differentially expressed circRNAs in MPE. KEGG analysis showed that there were 14 pathways related to the top 5 up-regulated circRNAs, including PI3K-Akt signaling pathway, Regulation of actin cytoskeleton, and Rap1 signaling pathway. PI3K-Akt signaling pathway plays an important role in tumorigenesis and regulates critical cellular functions including survival, proliferation and metabolism [49]. Regulation of actin cytoskeleton is an essential step in cell migration during activation of the epithelial-mesenchymal transition program, which is associated with metastatic properties of cancer cells [50]. And Rap1 signaling pathway has also been reported to be involved in tumorigenesis and tumor progression [51]. All these results strongly indicate that these top 5 up-regulated circRNAs may be of importance in the development and progression of MPE. At the same time, we also performed GO analysis, illustrating the biological process, cellular component, and molecular function of the target genes of these 5 up-regulated circRNAs. These GO results showed that target genes of candidate circRNAs were mainly enriched in GO terms of protein binding, plasma membrane, and extracellular exosome, further suggesting that circRNAs may be involved in MPE development by regulating miRNA and their target genes.

As a conclusion, we identified multiple circRNAs which are differentially expressed in LA-MPE for the first time, proposed that the circRNA-miRNA-mRNA axis may be the possible molecular regulatory mechanism of MPE, and hoped the study can provide a reference for the broad analysis of circRNAs in MPE. Nevertheless, some limitations must be considered in this study. First, despite the appreciable effort of extending samples, the sample size was still limited and future studies should endeavor to recruit larger number of pleural effusion samples to confirm our results. Second, though the predicted genes of circRNAs were all subjected to LA-MPE associated microarray data and intersection genes were selected for analysis, the exact molecular mechanisms of candidate circRNAs needed to be elaborated in follow-up studies. Third, though we detected the expression profile of circRNAs in cell pellets in pleural effusion, its expression in cell-free supernatant need to be explored. However, studies showed that exonic circRNAs is not stable in serum, with a half-life of less than 15s [24].

Disclosure of Potential Conflicts of Interest

The authors have no conflicts of interest to declare.

Acknowledgments

We thank all investigators who were involved in this study. We express our thanks to our colleagues Dong Wei, Chunxue Ran, and Congfeng Li.

The Corresponding author Ziliang Hou thanks his wife Shengjiao Wang for her dedication and love. The acceptance of this paper will be a gift for their unborn child.

Additional information

Funding

This work was supported by grants from the Science and Technology Commission Foundation of Tongzhou District, Beijing, China [KJ2017CX036-07].