Differential Expression Profiles and Functional Analysis of Circular RNAs of Sinonasal Inverted Papillomas

Keywords

Sinonasal Inverted Papilloma; Circular RNAs, Gene Expression Profile; Microarray; ceRNA.

Introduction

Sinonasal inverted papilloma (SNIP) is one of the most common benign tumors of sinonasal cavities and it accounts for 0.5-4% of sinonasal tumors and about 70% sinonasal papilloma [1]. It is originated from ectodemal schneiderian membrane of nasal cavity or paranasal sinuses [2, 3]. Compared with nasal polyps, the lesions of SNIP are firmer, bulkier and lack translucency. Its main clinicopathological features include the tendency to recur, easily malignant transformation and its aggressive growth which can cause local invasion and bone destruction [4, 5]. Nasal obstruction, rhinorrhoea, epistaxis, headache and anosmia are main symptoms of SNIP. The optimal treatment for SNIP is surgery which aimed at removal the tumors completely. The etiology and pathogenesis of SNIP are still unclear. Study about molecular mechanisms can help us better understand the development of SNIP and may provide a better treatment for SNIP. Noncoding RNAs (ncRNAs) are RNAs that do not encode proteins, which including small nuclear RNAs, transfer RNAs and ribosomal RNAs, as well as the more recently discovered long noncoding RNAs (lncRNAs), circular RNAs (circRNAs) and microRNAs (miRNAs), and they have critical regulatory roles in many biological functions [6, 7]. CircRNAs are a class of ncR-NAs which were first discovered in human in 1986 [8]. They are present in the cytoplasm of eukaryotic cells and are singlestranded circular molecules formed by pre-mRNA back-splicing through variable shearing processes [9-12]. CircRNAs are more stable compared with linear RNAs because they are covalently closed loop structures with no 5’ to 3’ polarity and polyadenylation which make them resistance to RNA exonucleases [13, 14]. Potential biological functions of circRNAs is still unclear, but increased researches reported that circRNAs can act as competing endogenous RNAs (ceRNAs) of miRNAs via miRNA response elements (MREs)[15] which can regulate transcription and gene expression [16]. Many studies have found that circRNA play an important role in many life processes including aging, tissue development, insulin secretion, and cell apoptosis [17]. Also, circRNAs are related to several diseases, such as cancer, autoimmune disease, cardiovascular diseases and neurological disorders [18, 19]. However, few research focus on the role of circRNAs in SNIP and the expression profiles of circRNAs in SNIP are still undiscovered. In this study, we explored the differential expression profiles between the SNIP and normal nasal mucosa tissue using RNA microarray. Then, quantitative real- time reverse transcription PCR (qRT-PCR) was used to validate the selected differentially expressed circRNAs.

Materials and Methods

Patients and Sample Collection

SNIP specimens and normal nasal mucosa tissues outside the tumor were obtained from eight patients undergoing endoscopic resection of inverted papilloma at Department of Otorhinolaryngology of Shandong Provincial Hospital Affiliated to Shandong First Medical University. The specimens were collected between June 2020 and May 2021. All SNIP cases were confirmed by postoperative pathology. The exclusion criteria were as follows: recurrent SNIP, combined with other neoplastic lesions of nasal cavity and paranasal sinus, accompanied by systemic malignant tumors or immune diseases. The basic clinical characteristics of eight patients with SNIP are shown in table 1. The tumor tissue and normal tissue of each patient were used as controls. Besides, six pairs of samples (SNIP_1-6) were used for microarray analysis, and three pairs of samples (SNIP_6-8, including one pair used for microarray) were used by qRT-PCR for microarray data verification. The research was improved by the institutional Ethics Committee of Shandong Provincial Hospital Affiliated to Shandong First Medical University.

RNA extraction and microarray analysis

Total RNA was isolated using RNeasy Total RNA Isolation Kit (Qiagen, GmBH, Germany)/ TRIzol reagent (Life technologies, Carlsbad, CA, US) according to the manufacturer’s instructions, and purified by using a RNeasy Mini Kit (Qiagen, GmBH, Germany). Total RNA was checked for a RIN number to inspect RNA integration by an Agilent Bioanalyzer 2100 (Agilent technologies, Santa Clara, CA, US). RNA samples of each group were then used to generate biotinylated cRNA targets for the Sino Human ceRNA array V3.0. The biotinylated cRNA targets were then hybridized with the slides. After hybridization, slides were scanned on the Agilent Microarray Scanner (Agilent technologies, Santa Clara, CA, US). Data were extracted with Feature Extraction software 10.7 (Agilent technologies, Santa Clara, CA, US). Raw data were normalized by Quantile algorithm, R package “limma”. The microarray experiments were performed by following the protocol of Agilent technologies Inc at Sinotech Genomics Corporation. Genes with a fold change of at least 2 were selected for further analysis. Heatmap plots were done by a R package “pheatmap” of the target genes.

Functional analysis by GO and KEGG pathway

Gene ontology (GO) (www.geneontology.org) is a kind of ontology which is widely used in the bioinformatics field. GO categories were divided into three parts according to its function: biological process (BP) and cellular component (CC), molecular function (MF). It was used to explore the functions of up-regulated and down-regulated genes. Kyoto encyclopedia of genes and genomes (KEGG) (www.kegg.jp) pathway analysis is a frequently used analysis of gene function and genome information database and it can help researchers to find out gene and expression information. KEGG enrichment analysis of dysregulated genes can find out differentially enriched pathways, which is helpful to identify the functions of these dysregulated genes in the biological regulation pathway. GO/pathway enrichment analysis were done use Fisher's exact test by a R package “clusterProfiler” of the target genes. GO categories/Pathway with Fisher's exact test P values < 0.05 were selected.

Quantitative real time-PCR

Quantitative real time-PCR (qRT-PCR) was carried out by using LightCycler 480 (Roche, Shanghai, China) system to quantitate the levels of three up-regulated circRNAs (has_circ_0025388, has_circ_0014213, has_circ_26421) and three down-regulated circRNAs(has_circ_0027511, has_circ_0081375, has_ circ_0013507). The qRT-PCR was conducted by the instructions of the SYBR® Green Premix Pro Taq HS qPCR Kit (AG, China). The primers sequences of circRNAs for qRT-PCR validation are shown in Table 2. The relative expression levels of circRNAs were presented using 2-ΔΔCt method.

---

---

---

---

Co-expression network

Co-expression network is to calculate the co-expression relationship between genes according to the dynamic change of gene expression signal value, obtain the expression regulation relationship and regulation direction between genes, and then construct the gene expression regulation network. In order to better understood the function of circRNAs and the relationship between circRNA and mRNA, we calculated the co-expression relationship between circRNA and mRNA, so as to construct the gene expression regulation network. By using the co-expression network, we can analyze the gene regulation ability and obtain the core regulatory genes that the samples change with the experiment. We construct the co-expression network between circRNAs and mRNAs using Cytoscape.

Competing endogenous RNA (ceRNA) network construction

Recent studies have shown that circRNA molecules are rich in microRNA binding sites, which can act as miRNA sponge in cells, and further relieve the inhibition of miRNA on its target genes, and increase the expression level of the target genes. This mechanism is called the competitive endogenous RNA (ceRNA) mechanism. The ceRNA network base on this hypothesis which reveals a new mechanism of RNA interaction that ceRNA (mRNA, lncRNA, circRNA) can regulate each other's expression by competitive binding miRNA [20, 21]. CeRNA analysis is based on gene expression values. Through regression model analysis and seed sequence matching, a regulatory network for sponge adsorption of microRNA is established to find the core ceRNA. Interactions of circRNAs with potential miRNA response elements (MREs) were assessed with Cytoscape software. Then we constructed the circRNA-miRNA-mRNA network of the two circRNAs (has_ circ_0026421 and has_circ_0027511) which verified in qRT-PCR.

Statistical analysis

The data were expressed as mean ± standard deviation (SD). The differences between the two groups of data were compared by t-test using SPSS 24.0 software. p<0.05 indicated statistically significance.

Results

Dysregulatated circRNA and mRNA

In order to explore the dysregulated circRNA and mRNA in SNIP, the microarray analysis was performed using SNIP samples and normal nasal mucosa samples. According to the results of circRNA microarray analysis, 4264 circRNAs were significantly dysregulated (fold change≥2, p<0.05) which contain 1567 upregulated and 2697 down-regulated (Fig. 1). Based on the result of mRNA microarray analysis, 4194 mRNAs were significantly dysregulated (fold change≥2, p<0.05) , including 1894 up-regulated and 2300 down-regulated (Fig. 2). Of which the top 10 upregulated and down-regulated circRNAs were shown in Table 3.

---

---

---

---

---

---

Function analysis of dysregulated circRNAs

GO analysis were carried in order to analyze which gene functions are related to differentially expressed circRNAs. The top 3 enriched GO terms of dysregulated circRNAs were “cellular process” “single-organism process” “biological regulation” in BP; “cell” “cell part” “organelle” in CC; “binding” “catalytic activity” “molecular function regulator” in MF (Fig. 3A). Similarly, KEGG pathway analysis was used to explore which pathways were related to dysregulated circRNAs. And the results of top 30 pathways between SNIP samples and normal nasal mucosa samples were shown in Fig.3B. The focal adhesion pathway, apelin signaling pathway and ECM-receptor interaction pathway were enriched in dysregulated circRNAs. These results suggested that these pathways related to regulation of cell growth, proliferation differentiation and degradation may contribute to the pathogenesis of SNIP.

Validation of circRNA expression

Six circRNAs including three up-regulated and three downregulated which were chosen from microarray data were validated using the qRT-PCR. Of which has_circ_0014213 and has_ circ_0025388 failed in qRT-PCR analysis because of nonspecific products. The expression trend of has_circ_0013507 was inconsistent. Among the remaining three circRNAs, there was no significant difference in the relative expression of has_circ_0081375 between normal and SNIP tissues. While the results of other two circRNAs (has_circ_0026421 and has_circ_0027511) were consistent with the microarray, and the relative expression differences were statistically significant (Fig. 4).

Co-expression network construction of circRNAs and mRNAs

According to the microarray data and calculation results of a Pearson’s correlation coefficient, numerous coordinately expressed genes’ networks were construction in SNIP tissue compared with the normal nasal mucosa. By analyzing the functional analysis (GO analysis and KEGG pathway analysis) of the top 20 circRNAs (ten most up-regulated and ten down-regulated), we selected two pathways of Rap 1 signaling pathway and cAMP signaling pathway. Based on these two pathways, we constructed a circRNA-mRNA co-expression network of the ten most up-regulated and down-regulated circRNAs (Fig. 5). The results shows that circRNAs work together with other circRNAs and mRNAs rather than act alone.

ceRNA network construction of circRNA-miRNA-mRNA

CircRNA molecules are rich in miRNA binding sites and play the role of miRNA sponge in cells, so as to relieve the inhibition of miRNA on its target genes and increase the expression level of target genes. This mechanism is called ceRNA mechanism. Considering microarray analysis and verification results, has_ circ_0026421 and has_circ_0027511 were expressed significantly differentially in SNIP. In order to explore the possible pathways of the verified circRNAs, based on the two pathways mentioned above (Rap 1 signaling pathway and cAMP signaling pathway), we constructed circRNA-miRNA-mRNA network of these two selected circRNAs which was shown in Fig. 6. The results showed that the top two MREs of hsa_circ_0027511 are hsa_miR-548v and hsa_miR-557, and the top two MREs of hsa_circ_0026421 are hsa_miR-6512-3p and hsa_miR-6722-3p.

---

---

---

---

---

---

---

---

Discussion

Sinonasal inverted papilloma is a benign epithelial tumor, which originates from Schneider mucosa of nasal cavity and paranasal sinus [22]. Although the exact etiology of SNIP remains unknown, it is a multifactorial disease that may be related to human papilloma virus (HPV) infection, chronic inflammation, allergy and occupational exposure [23]. At present, the diagnosis, treatment and follow-up of SNIP are mainly based on imaging, pathology and endoscopic technology. But there was still a high recurrence rate and malignant transformation rate after operation. With the development of molecular biology technology, the research on the mechanism of SNIP occurrence, development, proliferation and malignant transformation has made more in-depth progress at home and abroad. Li et al.[24] found that the overexpression of ErbB1 and ErbB2 in SNIP epithelial cells may play an important role in the pathogenesis and progress of SNIP. ErbB1 and erbB2 may be new targets for the treatment of SNIP. Yasukawa et al.[25] confirmed that kirsten rat sarcoma viral oncogene (KRAS) gene is a key factor in predicting the malignant transformation of SNIP in Japanese patients. Lin et al.[26] found that Stathmin was up-regulated in SNIP patients by immunohistochemical staining. Yamashita et al.[27] proposed that there was a significant correlation between serum squamous cell carcinoma antigen (SCCA) level and SNIP size, and the increase of SCCA was a significant predictor of tumor size. These studies may be helpful to explore the key genes affecting SNIP and provide new ideas and favorable evidence for the diagnosis and treatment of SNIP. However, most studies have not been used in clinic because of insufficient sample size, geographical limitations and inability to identify key factors. As is known to all, many human genome sequences do not encode any protein, and studies found that non-coding RNA accounts for almost 95% of the total RNA sequences transcribed from the human genome, including circular RNA[28]. At the time when the circRNA was first be discovered, it was considered to be a low-abundance RNA molecule which formed by mis-splicing of exons when transcription [29]. With the development of technology of bioinformatics and RNA sequencing, circRNAs are be found present widely in the eukaryotic cells and play an key role in regulating gene expression at the post-transcriptional level. So far, they have been reported to be related to regulating immunity, inflammations, and cell proliferation and they are thought to play essential roles in a number of diseases [30-33]. Based on the above research results, we considered that circRNAs may play an important role in SNIP. However, accurate expression profile of circRNAs and their possible functions in SNIP have not been reported so far. In this present study, we try to explore the expression profile of circRNA in SNIP and analyze its potential functions.

Firstly, we explored the expression profiles of circRNAs of SNIP tissue and normal nasal mucosa extracted from six patient with SNIP when undergoing endoscopic sinus surgery. As mentioned earlier, a total of 4264 circRNAs were significantly dysregulated (fold change≥2 , p<0.05) which contain 1567 up-regulated and 2697 down-regulated. Then, in order to confirm the accuracy of the microarray results, we selected 6 abnormally expressed circRNAs (including 3 up-regulated and 3 down-regulated), which were verified by qRT-PCR in three groups of tumor and adjacent normal tissue specimens. Three of them were confirmed to show the same trend of circRNAs in the qRT-PCR and microarray, in which, relative expression differences of two circRNAs were statistically significant, including one was up-regulated(has_circ_0026421) and the other was down-regulated (has_circ_0027511). Then, we constructed the co-expression network of circRNAmRNA and circRNA-miRNA-mRNA network to obtain the relationship between circRNAs and mRNAs. And we focused on the analysis of ten most up-regulated and down-regulated circRNAs, especially the two verified circRNAs (has_circ_0026421 and has_circ_0027511). Results showed that the top two MREs of hsa_circ_0027511 are hsa_miR-548v and hsa_miR-557, and the top two MREs of hsa_circ_0026421 are hsa_miR-6512-3p and hsa_miR-6722-3p. These finding can enable us to further understand the pathogenesis of SNIP.

The GO and KEGG analyses were performed to understand the biological functions of dysregulatated circRNA. The three most enriched GO terms were “cellular process” “single-organism process” and “biological regulation”. And KEGG pathways were enriched “focal adhesion pathway”, “apelin signaling pathway” and “ECM-receptor interaction pathway”. These pathways were related to regulation of cell growth, proliferation differentiation and degradation [34-36]. Thus, these current findings mean the dysregulated circRNAs might contribute to SNIP via these pathways and through these biological processes. Further researches will focus on the action sites and regulation modes of these pathways and biological processes.

In the future study, we will focus more on has_circ_0026421 and has_circ_0027511, and more experiments will be used to further evaluate the role potential role of has_circ_0026421 and has_ circ_0027511 functioned as ceRNA in SNIP.

Conclusion

In conclusion, this study demonstrated that circRNAs are aberrantly expressed in SNIP. Has_circ_0026421 was significantly up-regulated and has_cirC_0027511 were significantly down-regulated in SNIP tissues, which suggested that they may be novel biomarker for SNIP diagnosis and targeted therapy. In the future, more studies are needed to investigate the potential of circRNAs as biomarkers for SNIP diagnosis and targeted therapy. Funding

This work was supported by the Shandong Province Medical and Health Technology Development Plan Project (grant no.2019WS500), the Key Technology Research and Development Program of Shandong (Grant no. 2019GSF108257 and 2018GSF118192).