Cholesteatoma tissue and matched post-auricular normal skin tissue were sampled in each patient. A total of 7 patients aged 18–32 years who underwent unilateral middle ear cholesteatoma surgeries between June 2016 and November 2016 at the Department of Otorhinolaryngology, Peking Union Medical College Hospital were enrolled in this study. All participants provided written informed consent and were clinically and histologically confirmed to present with cholesteatoma. All specimens were stored at −80°C immediately after collection for later RNA extraction. Ethical approval for the study was obtained from the Ethics Committee of Peking Union Medical College Hospital (no. S-K292).

Total RNA was extracted from cholesteatoma and post-auricular skin tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer's instructions. RNA quantity and quality were assessed by NanoDrop ND-1000. RNA integrity was assessed by standard denaturing agarose gel electrophoresis or Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

Four pairs of cholesteatoma and post-auricular skin tissues were used for microarray assay to determine differentially expressed lncRNAs and mRNAs comparing cholesteatoma and post-auricular skin specimens. Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technologies, Englewood, CO, USA) with minor modifications. The hybridized arrays were washed, fixed, and scanned using the Agilent DNA Microarray Scanner (part no. G2505C). Agilent Feature Extraction software (version 11.0.1.1) was used to analyze acquired array images. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v12.1 software package (Agilent Technologies, Englewood, CO, USA). After quantile normalization of the raw data, lncRNAs and mRNAs in the 8 samples flagged as ‘Present’ or ‘Marginal’ were chosen for further data analysis. Differentially expressed lncRNAs and mRNAs with statistical significance between the two groups were identified through P-value/FDR filtering. Differentially expressed lncRNAs and mRNAs between the two samples were identified through Fold Change filtering. Pathway analysis and Gene ontology (GO) analysis were applied to determine the roles played by these differentially expressed mRNAs in these biological pathways or GO terms. Hierarchical clustering and combined analysis were performed by using in-house scripts.

GO provides a controlled vocabulary to describe gene function and relationships between these concepts in any organism (http://www.geneontology.org). GO covers three aspects: Biological process, cellular component, and molecular function. Fisher's exact test was applied to determine whether the overlap between the differentially expressed (DE) list and the GO annotation list is greater than that expected by chance. The -log₁₀ (P-value) was used to denote the significance of the GO term enrichment in the DE genes. FDR represents the false discovery rate. The lower the P-value is, the more significant the GO term is (a P-value <0.05 is recommended). Pathway analysis was performed to collect pathway clusters on the molecular interaction and reaction networks by mapping genes to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (http://www.genome.jp/kegg/). The -log₁₀ (P-value) denotes the significance of the pathway correlations. The lower the P-value is, the more significant the correlation is (a P-value <0.05 is recommended).

The selected lncRNAs and primers used for qRT-PCR were designed using Primer 5.0 and synthesized by Generay Biotech (Shanghai, China). β-actin was used as an internal control for all samples. Primer sequences were as follows: β-actin forward, 5′-GTGGCCGAGGACTTTGATTG-3′ and reverse, 5′-CCTGTAACAACGCATCTCATATT-3′; ENST00000415386 forward, 5′-TGGAGTAGGCACAGACGGAA-3′ and reverse, 5′-GACTGTGGTACAATCGCTTCG-3′; ENST00000420253 forward, 5′-CGATGTCAACCCTCGAACCT-3′ and reverse, 5′-TCAGATGTGGCCCACTGTCT-3′; NR_024468 forward, 5′-GCTTACATTTTCTCTGCCATCTC-3′ and reverse, 5′-CCACTTTCCTTTCTCTCTCTTCT-3′; T044224 forward, 5′-TGCGAAAACACTGCGATTG-3′ and reverse, 5′-GGTTGGACCAGACCAGATGA-3′; T347175 forward, 5′-GAAAGAGGCACTGCTGTTGA-3′ and reverse, 5′-GGCTGCTCCCAGAATAGATAG-3′; uc001kfc.1 forward, 5′-TTCCCAGAAGGCGTAGGTT-3′ and reverse, 5′-GACAGAGTCTCCCTCTATCATCC-3′; qRT-PCR was performed by using ViiA 7 Real-time PCR System (Applied Biosystems, San Jose, CA, USA) with a SYBR expression assay system (Takara Biotechnology, Co., Ltd., Dalian, China). The PCR reaction conditions were as follows: An initial denaturation at 95°C for 10 min, followed by 40 PCR cycles at 95°C for 10 sec and 60°C for 60 sec. Then annealing and extension at 95°C for 10 sec, 60°C for 60 sec and 95°C for 15 sec. Each sample was assayed in triplicates. The 2^−ΔΔCt method was used to determine fold-change in gene expression in the cholesteatoma samples relative to the normal skin samples. For statistical analysis, we used unpaired t-test to compare the expression of lncRNAs between cholesteatoma and normal skin samples. P<0.05 was considered to indicate a statistically significant difference.

All potential miRNA response elements were searched based on the sequences of lncRNAs and mRNAs. LncRNA/miRNA/mRNA interactions were predicted by the overlapping of the miRNA seed sequence binding site both on the chosen dysregulated lncRNAs and the significantly dysregulated mRNA. miRBase V19 was selected to interact with the lncRNAs and mRNAs (http://www.mirbase.org/). LncRNA/miRNA interactions were predicted by miRcode (http://www.mircode.org/). miRNA/mRNA interactions were predicted by miRanda (www.microrna.org/) and TargetScan (http://www.targetscan.org/).

A microarray analysis was performed to profile differences in lncRNA and mRNA expression between 4 pairs of cholesteatoma and matched normal skin samples. According to microarray expression profiling data, 11,815 lncRNAs and 7,692 mRNAs were detected. The lncRNAs were carefully collected from the most authoritative databases such as RefSeq, UCSC Known Genes, Gencode, and other related sources. According to their relation with protein-coding genes, all deregulated lncRNAs in cholesteatoma were classified into six categories: Bidirectional (4.46%), exon sense-overlapping (1.16%), intron sense-overlapping (3.47%), natural antisense (8.75%), intronic antisense (12.71%) and intergenic (69.47%) (Fig. 1). Hierarchical clustering revealed that lncRNA and mRNA expression patterns between cholesteatoma and matched normal skin tissues were distinguishable (Fig. 2A and B). Volcano plots were used for visualization and assessment of the variation (or reproducibility) of lncRNA and mRNA expression between cholesteatoma and matched normal skin tissues (Fig. 2C and D).

Furthermore, by setting a threshold for differential expression at changes ≥2.0-fold, we identified 787 lncRNAs and 591 mRNAs that were differentially expressed (P<0.05) between cholesteatoma and matched normal skin tissues (Fig. 2C and D). Among them, in cholesteatoma samples, 181 lncRNAs and 155 mRNAs were upregulated (fold change ≥2.0, P<0.05) and 606 lncRNAs and 436 mRNAs were downregulated (fold change ≥2.0, P<0.05) compared with normal skin samples (Fig. 2C and D). The microarray profile and RNAseq data sets have been deposited into Gene Expression Omnibus (GEO) with accession number GSE102673 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE102673).

GO analysis was performed to determine the gene and gene product enrichment in categories including biological process (BP), molecular function (MF), and cellular component (CC). LncRNAs can regulate the neighboring and overlapping coding gene expression. Thus, GO enrichment analysis of differentially expressed mRNAs may partially display the role of differentially regulated lncRNAs. In the present study, the target genes of deregulated mRNAs were analyzed and compared with those in cholesteatoma. In the GO analysis, by setting P<0.05, the upregulated and downregulated mRNAs were analyzed separately and the top-10 enriched GO terms of BP, CC, MF were listed (Fig. 3). In BP analysis, regulation of phosphatidylinositol 3-kinase signaling (GO:0014066), phosphatidylinositol 3-kinase signaling (GO:0014065), and positive regulation of phosphatidylinositol 3-kinase signaling (GO:0014068) belonged to top-10 most enriched processes associated with the upregulated mRNAs. Moreover, phosphatidylinositol 3-kinase binding (GO:0043548) was the most significantly enriched function of the upregulated mRNAs in MF analysis (Fig. 3C). It is notable, since phosphatidylinositol 3-kinase (PI3K) has been reported to play a crucial role in the formation of cholesteatoma (24). Furthermore, growth factor activity (GO:0008083), one of the most significantly enriched processes associated with the downregulated mRNAs in MF, is also involved in cholesteatoma formation (25).

By mapping genes to KEGG pathways, we performed pathway analysis (http://www.genome.jp/kegg/pathway.html). KEGG pathway enrichment analysis for differentially expressed mRNAs is devised to comprehend pathways and molecular interactions related to genes. Pathway analysis indicated that 10 enriched pathways corresponded to downregulated mRNAs and 7 pathways corresponded to upregulated mRNAs (P<0.05, Fig. 4). The pathways enriched with upregulated lncRNAs were involved in a category ‘bacterial invasion of epithelial cells (hsa05100)’ that is involved in the molecular mechanisms of cholesteatoma (26). Furthermore, ‘viral myocarditis (hsa05416)’ and ‘Hepatitis B (hsa05161)’ were also related to the molecular biology of cholesteatoma (6).

Using qRT-PCR, 3 upregulated and 3 downregulated lncRNAs with fold-changes >2.0 were randomly selected to verify the microarray data in 7 pairs of samples (4 pairs of original tissues and 3 other pairs of cholesteatoma and normal skin tissues). The qRT-PCR results and microarray data were consistent (Fig. 5), demonstrating that the microarray expression results are highly reliable.

CeRNAs are involved in a regulatory mechanism between non-coding RNA and coding RNA based on shared MREs (12). According to the ceRNA hypothesis, ceRNA members can compete for the same MREs to regulate each other. To explore whether lncRNAs have ceRNA potential in the pathogenesis of cholesteatoma, in the present study, we constructed a ceRNA network in cholesteatoma based on the microarray data. We selected 5 significantly differentially expressed lncRNAs (fold change >2.0, P<0.05, Table I), which shared common binding MREs with each other, to constitute an lncRNA/miRNA/mRNA ceRNA network (available upon request). The network was found to be composed of 20 lncRNA nodes, 399 miRNA nodes, and 137 mRNA nodes. In the network, we observed that lncRNA-uc001kfc.1 interacted with miR-21-3p (Figs. 6 and 7A), a microRNA belonging to the miR-21 family, which promotes the formation and invasion of cholesteatoma (7,8). (To make the lncRNA-uc001kfc.1-mediated ceRNA network easier to identify, we separated it from the original ceRNA network). The 2D structure of the miR-21-3p on uc001kfc.1 indicated that the binding site is 8 mer, perfectly matching positions 2 to 8 of the mature miR-21-3p (Fig. 7A).

Furthermore, to better explore the ceRNA potential of lncRNAs in cholesteatoma, we performed GO and KEGG pathway analysis for the lncRNA/miRNA/mRNA network. In GO analysis, gene product enrichment in BP, MF, and CC were analyzed (P<0.05) and the top-10 enriched GO terms of the three aspects are listed in Fig. 7B-D. GO analysis revealed that protein kinase activator activity (GO:0030295) as well as protein kinase A binding (GO:0051018) and protein kinase A regulatory subunit binding (GO:0034237) were the most significant molecular functions displayed in the ceRNA network. Protein kinase can regulate protein phosphorylation, which is involved in many aspects of cell biology and contributes to cholesteatoma formation (27). In KEGG pathway analysis, as depicted in Fig. 7E, the cAMP signaling pathway, the JAK/STAT signaling pathway, and the mTOR signaling pathway were related to cell migration, cell cycle, and survival (http://www.genome.jp/kegg/). Moreover, JAK/STAT signaling was also involved in the molecular mechanisms underlying cholesteatoma (28).