Clinical biomarkers of pulmonary carcinoid tumors in never smokers via profiling miRNA and target mRNA

Background miRNAs play key regulatory roles in cellular pathological processes. We aimed to identify clinically meaningful biomarkers in pulmonary carcinoid tumors (PCTs), a member of neuroendocrine neoplasms, via profiling miRNAs and mRNAs. Results From the total of 1145 miRNAs, we obtained 16 and 17 miRNAs that showed positive and negative fold changes (FCs, tumors vs. normal tissues) in the top 1% differentially expressed miRNAs, respectively. We uncovered the target genes that were predicted by at least two prediction tools and overlapped by at least one-half of the top miRNAs, which yielded 44 genes (FC<-2) and 56 genes (FC>2), respectively. Higher expressions of CREB5, PTPRB and COL4A3 predicted favorable disease free survival (Hazard ratio: 0.03, 0.19 and 0.36; P value: 0.03, 0.03 and 0.08). Additionally, 79 mutated genes have been found in nine PCTs where TP53 was the only repeated mutation. Conclusion We identified that the expressions of three genes have clinical implications in PCTs. The biological functions of these biomarkers warrant further studies.


Background
Pulmonary carcinoid tumors (PCTs) are a member of lung neuroendocrine neoplasm family. The incidence of PCTs is relatively low and comprises approximately 1-2% of lung neoplasms [1][2][3]. PCTs occur frequently in never-smokers, and its molecular etiology is still unclear. PCTs are subdivided into typical carcinoids (TCs) and atypical carcinoids (ACs), where TCs lack necrosis and have < 2 mitoses/2 mm 2 , and ACs have 2-10 mitoses/2 mm 2 and/or necrosis [3]. Seventy to ninety percent of PCTs are identified as TCs and 10-30% as ACs [4]. The overall survival of cases undergoing radical resection are relatively satisfactory: 92-100% for TCs and 61-88% for ACs [3,5]. However, surgical resection is still the only effective treatment, and unresectable tumors are very difficult to treat due to their insensitivity for both radiation and chemotherapy. Additionally, recurrence or metastasis can occur even decades after primary resection [6]. miRNAs are genomically encoded small non-coding RNAs that regulate genetic information by controlling stability or translation of mRNAs. Thus far, the study on exploring miRNA biomarkers and their target genes for disease monitoring or as a treatment option is state-ofthe-art, timely and highly promising [7,8]. We believe the profiling of miRNAs and screening of their target mRNAs are important to elucidate the tumorigenesis mechanisms of PCTs and indicate the potential directions of targeted therapy. In the study, we investigated the differential expression of miRNAs and target mRNAs. Additionally, we analyzed the clinical implications of the expression of these miRNAs and mRNAs in PCTs.

Results
Fifty-one individuals who had never smoked tobacco products and developed PCTs between 1997 and 2008 were studied. All study protocols were approved by Mayo Clinic's Institutional Review Board. The clinical and demographical features of these cases are shown in Table 1.

Sources of variation (SOV) and principal components analysis (PCA) of miRNAs expression
As indicated in the Additional file 1: Figures S1A and B, tumor status (tumor vs. normal) was the most effective factor with regard to SOV of miRNAs expression. Plate effects that were resulting from the split of the testing into manageable rounds, i.e., batch effects within each sample type, had the second highest significant F ratio. Hence, we used the plate effects as the covariance in the ANOVA model for differential analysis between tumor and normal.
As indicated in Additional file 1: Figures S1C and D, PCA indicates that miRNA expressions between tumors vs. normal tissues of formalin-fixed paraffin-embedded tissues (FFPE) were better separated than fresh-frozen (FF) samples.
The profiling of miRNAs in the top 1% of fold changes (FCs) We identified the median FC, Bonferroni P-value and False Discovery Rate (FDR)-Q-value for all of the 1145 miRNAs. We sorted miRNAs by FC in FFPE and FF sample sets, respectively. Thereafter, we selected the top 12 (12/1145, i.e., 1%) in FF and FFPE, respectively. Table 2 lists the 16 miR-NAs that showed a positive FC. Among the 16 miRNAs, eight miRNAs are the top 12 in both of the sample sets, and the other eight miRNAs are the top 12 in either of the sample sets. Table 3 lists the top 17 miRNAs that showed a negative FC. Among the 17 miRNAs, seven miRNAs are the top 12 in both of the sample sets, and the other ten miRNAs are the top 12 in either of the sample sets. Table 2 indicates that FCs ranged from 10.86 to 5.51 and from 14.04 to 6.19 in FF and FFPE sample sets, respectively. Table 3 indicates FCs ranged from -5.28 to -2.08 and from -3.47 to -1.75 in FF and FFPE sample sets, respectively. The absolute value of positive FCs were remarkably greater compared to the absolute values of negative FCs, suggesting PCTs are prone to have inhibited tumor suppressor genes via high expression of miRNAs.
The profiling of target mRNAs with FCs>2 or <-2 To screen the target mRNAs of the miRNAs as in the top 1% of FCs, we used prediction tools to predict the target mRNAs of the miRNAs listed in Tables 2 and 3, respectively. Thereafter, we uncovered the target genes which were predicted by at least 2 prediction tools and overlapped by at least one-half of the miRNAs in Table 2 (8 of  16 miRNAs) and Table 3 (9 of 17 miRNAs), respectively.
The expressions of 680 and 792 target genes, which were derived from the miRNAs in Tables 2 and 3, were finally confirmed in the DASL (cDNA-mediated Annealing Selection Extension and Ligation) assay, respectively. We listed all 44 genes with a median FC<-2 (expressions in tumors are at least two-fold lower than in normal tissues, Table 4) and 56 genes with a FC>2 (expressions in tumors are at least two-fold higher than in normal tissues, Table 5), respectively.

The expressions of miRNAs and mRNAs in tumors and disease free survival (DFS)
In order to analyze the association between expression levels of miRNAs, mRNA, and clinical outcomes, we combined FF and FFPE together and adjusted sample type as a covariate to calculate the hazard ratio (HR) due to the small numbers of each sample set. Thus, there were a total of 47 cases in the miRNA group and 48 cases in the mRNA group (one patient more in the mRNA group compared to the miRNA group) who had available follow up information. There were five cases with status=1 (death, recurrence or progression). The median follow up time was 6.4 years. All the cases underwent surgery, and one case underwent radiation therapy following surgery. There was no correlation between DFS and clinical-demographical characteristics, i.e., gender, age and stage.

miRNAs and DFS
We used the expression levels of the miRNAs, which are listed in Table 2 and as the predictors, respectively. However, there was no miRNA with significant predictive power for DFS.

mRNAs and DFS
We analyzed the predicted mRNAs in tumors, which are listed in Tables 4 and 5. Among the genes with a FC<-2, two genes, i.e., CREB5 and PTPRB had a significant HR, and one gene, COL4A3, had a marginally significant HR, suggesting high expressions of these genes remarkably predicted better DFS (Table 6). Additionally, no genes with a FC>2 had significant predictive power for DFS.
Validation of differential expressions of CREB5, PTPRB and COL4A3 by using RNA-sequencing (RNA-Seq) To validate the differential expression of CREB5, PTPRB and COL4A3, we used RNA-seq to detect the expression of the abovementioned genes in 10 pairs of tumor/ normal fresh frozen tissues. Additional file 1: Table S1 shows the Bonferroni-P-values, FDRs (Q-values) and FCs of the abovementioned genes by using DASL (FF sample set) or RNA-seq. All FCs of the mRNAs in DASL can be validated by RNA-seq (Additional file 1: Table S1).
We also evaluated the expression data of CREB5, PTPRB and COL4A3 from Oncomine [9], indicating that CREB5, PTPRB and COL4A3 were significantly lower in carcinoid than normal lung tissues (Additional file 1: Figure S1), which was accordant to the results of the DASL and RNAseq datasets.

Discussion
The molecular mechanism of the tumorigenesis of PCTs needs to be further investigated. Furthermore, there are very few studies focused on PCTs and miRNAs [18], which hold key roles in the regulation of a variety of cellular processes. It is well known that dysregulations of miRNAs are important to keep the malignant phenotype of many tumors. The distorted and unique expression profile of miRNAs in different types and subsets of tumors make miRNAs an attractive source of sensitive biomarkers [19]. Importantly, miRNAs not only regulate the expression of oncogenes and tumor suppressors but also play the roles as oncogenes and tumor suppressors. It is very important to screen the expression of these miRNAs in tumor and normal tissues by using miRNA profiling. In our study, we compared the differential expression of 1145 miRNAs between carcinoid tissues vs. matched normal lung tissues, by using miRNA microarray, a high throughput approach. Recently, Kolbert et al. [20] compared the results of miRNA expression levels by using several platforms, including Illumina miRNA arrays and Illumina Next Generation Sequencing, and verified these results with those of quantitative PCR data. It was demonstrated that the within-platform reproducibility for each method was consistently high and detection of miRNA transcripts was similar across multiple platforms. As a result, we believe it is reliable to test the levels of miRNA expression by using Illumnia arrays.
A major shortcoming of gene expression profiling studies on neuroendocrine lung tumors is the absence of a proper reference tissue because pulmonary neuroendocrine cells are only scarcely present in the lung. Therefore, normal adjacent lung tissue is used as a reference.
We obtained the top 16 miRNAs that showed positive FCs and 17 miRNAs with negative FCs, confirmed by either FF or FFPE sample sets. Although we did not identify the top miRNAs with significant predictive power for DFS probably due to the small sample size, we believe that the biological functions of the other top miR-NAs in PCTs still warrant further investigations.  By bioinformatics tools, we predicted the target mRNAs, which were derived from the miRNAs with top positive FCs and negative FCs, respectively. Finally, we obtained 44 genes with FC<-2 and 56 genes with FC>2, respectively.
The COX proportional hazard model demonstrated that doubled expressions of three genes, i.e., CREB5, PTPRB and COL4A3, was associated with 63% to 97% decreased risk of death, recurrence or progression, suggesting their potential role as tumor suppressors ( Table 6). The statistical power of the COX model of three genes were 0.990 (CREB5), 0.611 (PTPRB) and 0.635 (COL4A3), respectively. Therefore, we will confirm the clinical implications of these genes in a large cohort of cases. PTPRB has been reported to be a drug addiction-associated gene [21], and CREB5 to be upregulated in the blood of cases with cluster headaches when compared to controls [22]. High COL4A3 expression was found to correlate with poor prognosis after combined cisplatin-gemcitabine chemotherapy in non-small cell lung cancer, showing the discrepancy of biological functions of this gene in different tumors [23]. We believe the above mentioned genes, which have remarkable differential expression and significant clinical implications in PCTs, to be potential therapeutic targets. However, the biological functions and clinical utility of these genes warrant further in vitro and in vivo investigations.
Recently, a study of 48 small intestine neuroendocrine tumors by parallel exome sequencing indicated that somatic single nucleotide variants affected a preponderance of cancer genes, including FGFR2, MEN1, HOOK3, EZH2, MLF1, CARD11, VHL, NONO and SMAD1 [24]. Our study indicates TP53 was the only repeated mutation (2/9), and each of the other 78 genes including HOOK3 were only mutated in one case, probably due to the disparities in different organs. Indeed, the p53 mutation has been reported in 20% of typical PCTs [25], which is in accordance to our results (2/9 cases). Intriguingly, some published studies [10][11][12][13][14][15][16][17] have proved that the mutation of p53 may regulate the expression of the miRNAs in a variety of tumors. We postulate that the p53 mutation may also affect the expression of the miR-NAs in PCTs. The underlying associations between the p53 mutation and miRNAs expression in PCTs warrant further studies.

Conclusion
We identified that the expressions of three genes have clinical implications in PCTs. The biological functions of these biomarkers warrant further studies.

miRNA expression assay
Four samples were excluded due to low quality. A group of 47 carcinoids tumors were analyzed; each with paired tumor-adjacent normal tissue samples (at least 5 cm away from the primary tumors). There are 24 pairs of tumor/ normal samples in FF tissues, and 23 pairs of tumor/ normal samples in FFPE. The clinical and demographical features of these cases are shown in Table 1. The expression status of 1145 miRNAs was generated by Illumina miRNA Arrays (San Diego, CA, U.S.A.), as described before [20]. We deposited the miRNA profiling data to Gene Expression Omnibus (GEO, accession number: GSE58600).
To explore the technical and biological factors that may affect miRNA expression levels, we performed multiple factor ANOVA analysis to identify SOV. To visualize the miRNA expression in tumor and normal tissues, we projected the numeric data by using scatter plot, i.e., PCA.

mRNA expression assay
A total of 51 carcinoid tumors were analyzed by the whole-genome DASL Assay (Illumina Beadchips Array, San Diego, CA, U.S.A.) as described before [26], each with paired tumor-adjacent normal tissue. There are 25 pairs of tumor/normal in FF tissues, and 26 pairs of tumor/normal samples in FFPE tissues. Additionally, nine pairs of tumor/ normal samples in FF tissues were analyzed by whole transcriptome sequencing (RNA-seq) using the TruSeq protocol (Illumina, San Diego, CA, U.S.A.). The clinical and demographical features of these cases are shown in Table 1. Finally, the differential expressions of mRNAs were validated by Oncomine (a public cancer microarray database; https://www.oncomine.org/).

Whole exome sequencing (WES) for gene mutation identification
Genomic DNA was purified from 9 pairs of carcinoid tumors and matched normal tissues. The nine cases were selected from the 51 cases that were analyzed for mRNA expression profiling. The clinical and demographical features of these cases are shown in Table 1

Survival analysis
Follow-up was conducted through detailed medical records data abstraction and self-administered questionnaires, starting within six months post diagnosis and then annually, thereafter. For deceased patients, the follow-up questionnaire was sent to the next-of-kin to acquire the death information, e.g., death date. We used disease free survival (DFS) as progression-free or recurrence-free survival from the operation. Recurrence, progression or death were counted as events. We evaluated the statistical power for survival analysis, considering numbers of events and sample size and Hazard ratio (HR).

Data analysis
The statistical analyses were performed using the analysis of variance (ANOVA) to analyze differential expression of miRNA and mRNA between tumors and normal tissues. Fold changes (FCs) of miRNA and mRNA differential expression were examined between tumors vs. normal tissues using stringent thresholds controlling for false discovery rates. Hence, a positive FC meant the expression level in