Introduction

Allergic rhinitis (AR) is a ubiquitous, lifelong illness that affects a large proportion of population and can be triggered by many allergens. This illness poses a heavy medical burden and can result in a poor quality of life.¹ Epidemiological studies have revealed that the prevalence of AR has increased progressively in recent decades with the rapid economic development and the changes in the lifestyles and dietary habits of people, especially children and teenagers.^2,3 The differentiation imbalance between T helper type 2 (Th2) and T helper type 1 (Th1) cells, both of which are CD4⁺ helper T lymphocytes, contributes to the pathogenesis of AR. The nasal mucosal inflammation of AR begins with allergen exposure and continues with the infiltration of B cells, T cells, basophils, eosinophils, and mast cells.⁴ However, the pathophysiology of AR remains ambiguous, and the currently available treatments often do not optimally control its symptoms.⁵ Therefore, an in-depth investigation of the new etiologies and treatment strategies for AR must be conducted.

The circular RNAs (circRNAs) are covalently linked RNA molecules generated via “back-splicing.”⁶ These molecules show an excellent stability, evolutionary conservation, and tissue- specific expressions. Given these properties, many studies have examined the potential role of circRNAs as disease biomarkers or therapeutic targets.⁷ Recently, numerous circRNAs have been identified as consequences of the rapid developments in high-throughput sequencing and bioinformatics. Given their function as miRNA sponges that regulate downstream mRNA expression and encode proteins, circRNAs have been examined in a growing number of studies by various forms of epigenetic modification.^8,9

An increasing amount of evidence reveals that circRNAs participate in the regulation of various immunocytes and immune responses. For instance, 597 circRNAs in CD4⁺ T cells were altered between asthmatic patients and controls.¹⁰ Given the many similarities in the inflammation profiles of AR and asthma, an aberrant circRNAs expression may also be found among AR patients. circRNAs may influence the polarization and differentiation of macrophage.¹¹ Some studies also highlighted the diverse functions of circRNAs in other immune cells, such as CD8⁺ T lymphocytes and humoral immune B cells.^12,13 Overall, the roles of circRNAs in immune-mediated diseases have been elucidated based on RNA-seq analyses and large-scale microarrays, including systemic lupus erythematosus, rheumatoid arthritis, and multiple sclerosis.^14–16 Nonetheless, the molecular mechanisms by which circRNAs regulate AR remain unclear, thereby calling for further investigations into this field.

In this study, the circRNA expression profiles in the nasal mucosa of a murine model were identified via RNA-Seq technology. The qRT-PCR was also performed to confirm the altered circRNAs. Several bioinformatic tools were used to explore the underlying functions of these circRNAs and to establish a circRNAs-miRNAs interaction network. Results show that altered circRNAs may be associated with the pathogenesis of AR and provide novel insights into this disease.

Materials and methods

Ethics statement

All animals used in this study were treated according to the animal experimental procedures permitted by the Animal Care and Use Committee of Shanghai Children’s Medical Center. Great efforts were made to reduce the number of experimental animals and alleviate their suffering. All animals were housed under standard specific pathogen-free laboratory conditions, and had free access to water and food. All mice were sacrificed after anesthetized deeply at the end of experiment.

Animal model

For the study, 6-week-old, healthy male BALB/c mice were obtained from SLAC Laboratory Animals center (Shanghai, China). These mice were randomly divided into the AR group (n = 4) and control group (n = 4). The BALB/c mice in the AR group were primary sensitized by intraperitoneal injection with a 200-μl mixture containing 100 μg OVA (Sigma- Aldrich, St. Louis, MO, USA) and 4 mg aluminum hydroxide (Sigma-Aldrich) on days 1, 8, and 15 of the experiment. From days 22 to 28, the AR group was challenged with intranasal drops of 20 μl OVA at a concentration of 40 mg/ml daily, whereas the mice in the control group were injected and challenged intranasally with a normal saline (NS) solution at the same schedule.

Measurement of induced allergic symptoms and sample collection

The symptom scores were evaluated after the final intranasal OVA or NS challenge on day 28. The sneezing, rhinorrhea, and nasal rubbing frequencies in each group were counted by blinded observers. Symptom scores were recorded based on the degree of severity and calculated as follows. Sneezing: the frequencies <3 were scored 1; 4–10 were scored 2; >11 were scored 3. Nasal rubbing: slight and occasional nasal rubbings were scored 1; repeated nasal rubbings were scored 2; rubbings from nose to face were scored 3. Rhinorrhea: observable watery discharge within the nasal cavity were scored 1; watery discharge spilling out of the anterior naris were scored 2; the face covered with abundant watery discharge were scored 3. A successful mouse model was determined by accumulative symptom scores >5. A day after the last challenge with OVA/NS, blood samples were collected from the retro-orbital venous plexus of the mice after deep anesthesia. Serum and plasma were obtained following a rapid centrifugation, and used for further enzyme-linked immunosorbent assay (ELISA) to measure the total immunoglobulin E (IgE) levels. After sacrifice, nasal mucosa was collected from the two groups of mice by using a small curette, and the collected specimens were immersed in liquid nitrogen and stored for further RNA-seq and qRT-PCR.

RNA isolation and sequencing

The total RNA from the nasal mucosa of the AR and control groups was isolated with a TRIzol reagent (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA) following the protocol of the vendor. The integrity and concentration of the extracted total RNA were then assessed by using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) and Qubit 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA), respectively. RNA samples with an RNA integrity number (RIN) value ≥7 and 28S:18S ratio ≥1.5 were used for subsequent processing. The NEBNext^® Ultra^TM Directional RNA Library Prep Kit for Illumina^® (NEB, USA) was utilized to generate a sequencing library. Library sequencing was performed on an Illumina HiSeq 3000 system with a pared-end program. The expression level of individual circRNAs was determined by calculating the reads per million mapped reads.

Fold change (≥2.0) and P value (<0.05) were used to identify the aberrant expressed circRNAs.

qRT-PCR

cDNA samples were synthesized from the total RNA via reverse transcription by using SuperScript^TM III Reverse Transcriptase (Invitrogen). The circRNA expression level was detected by using a SYBR Green PCR Master Mix (Applied Biosystems, Foster City, CA, USA). The qPCR reaction was achieved on an ABI 7500 real-time PCR system (Applied Biosystems): one 10-min cycle at 95 °C denaturation, followed by 40 cycles of 15 s at 95 °C and 60 s at 60 °C. The primers used to detect randomly chosen circRNAs to verify the results are shown in supplementary Table S1. The relative circRNA expression levels were calculated by using the 2^ΔΔCt method, and GADPH was used as internal control.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analysis

A GO analysis of the parental genes of the altered circRNAs was performed with a specific focus on biological processes, cellular components, and molecular functions. Pathway analysis was performed to manifest the diverse biological pathways of differential parental genes according to the KEGG database.

Prediction of circRNA-miRNA interaction

To construct a circRNA–miRNA network, circRNAs-targeted miRNAs were investigated by using TargetScan and miRanda-3.3 based on seed-match sequences. According to the association between the differentially expressed circRNAs and related miRNAs, circRNA–miRNA interaction network was generated and illustrated using Cytoscape 3.7.2.

Statistical analysis

All data were expressed as mean ± standard deviation, and were calculated in three independent experiments. All statistical data were analyzed by using SPSS version 22.0 and GraphPad Prism 8.0.1. The differential expression of circRNAs in the AR and control groups was tested by performing Student’s t tests. P < 0.05 was considered statistically significant.

Results

Murine model of AR evaluated by symptom scores and OVA-specific IgE

In AR studies, the symptom scores are the gold standard of a successful AR model. As shown in supplementary Figure S1A, nasal symptom scores were significantly higher in the AR group than in the healthy control group (2.52 ± 0.50 vs. 9.51 ± 1.03; P < 0.001). Since AR is caused by an IgE-mediated inflammatory, it is very important to show serum IgE levels. In Figure S1B, the concentration of OVA-specific IgE was significantly elevated in AR group compared to the normal group (7.78 ± 2.07 vs. 25.05 ± 2.98; P < 0.001). We came to the conclusion that OVA-sensitized AR model had been successfully established.

Figure S1 Establishment of OVA-induced AR murine model. (A) Total frequencies of sneezing, rubbing, and rhinorrhea were significantly increased in AR group compared with the healthy control group (P < 0.001). (B) OVA specific IgE in serum was markedly elevated in the AR group than in the control group (P < 0.001). ***indicates P < 0.001.

General characteristics of circRNAs

Secondary sequencing was performed to identify the circRNA expression profiles in the AR and control groups. A preliminary analysis was then performed according to the sequencing results. Box plots can be used to quickly visualize the distribution of a dataset for the circRNA profiles. After normalization, a similar distribution was observed among the eight samples (Figure 1A). These overlapped circRNAs were located throughout the whole genomic regions (Figure 1B), and the sequence length of circRNAs mostly ranged between 200 and 1600 nucleotides (Figure 1C). Among the screened-out 1103 circRNA candidates, 918 were already included in the circBase database and 185 were first identified. These circRNAs were distributed in each chromosome, including the X and Y chromosomes. Most of the identified circRNAs were transcribed from chromosomes 2 (8.25%), 5 (7.61%), and 9 (6.80%), whereas the smallest amount of circRNAs were transcribed from chromosome Y (0.18%). Most circRNAs were exonic and only few of them were intragenic, intronic, or antisense (Figure 1D). The raw RNA-seq data have been submitted to the National Center for Biotechnology Information Sequence Read Archive (SRA) under the SRA accession number PRJ609719 (https://www.ncbi.nlm.nih.gov/sra/PRJNA609719).

Figure 1 General characteristics of circRNAs identified in the AR group and the control. (A) Normalized RNA-seq data were presented as a box plot to visualize the distribution of circRNAs profiles. C1–C4 represented the samples in control group. A1–A4 represented the samples in AR group; (B) The distribution of total identified circRNAs on different chromosomes; (C) The length distribution of total identified circRNAs; (D) Classification of dysregulated circRNAs.

Differentially expressed circRNAs in AR murine models and controls

Hierarchical cluster analysis was performed to screen out expression patterns of circRNAs among the AR and control samples. A visual heat map showed that the expression patterns of circRNAs were distinguishable between the two groups (Figure 2A). The results of the RNA-seq analysis showed that 86 distinctly expressed circRNAs were filtered out, among which 51 were upregulated and 35 were downregulated (fold change ≥ 2.0, P < 0.05). The variation of circRNAs between the two groups were visualized via volcano plot filtering (Figure 2B). Table 1 summarizes the top 10 upregulated and top 10 downregulated circRNAs, among which chr5(+):17811180_17820600, chr9(+):111188761_111205847, and chr3(-):116778618_116782643 are novel circRNAs identified in this study. To verify the RNA-seq data, 10 differentially expressed circRNAs, including five upregulated and five downregulated circRNAs (nmu_circ_0005368, nmu_circ_0004409, nmu_circ_0001097, nmu_circ_0000950, nmu_circ_0013306, nmu_circ_0002130, nmu_circ_0004141, nmu_circ_0009939, nmu_circ_0014834, and nmu_circ_0000865), were randomly selected. The expression levels of these circRNAs were validated by performing qRT-PCR in the nasal mucosa of the murine AR model and the control group. The expression levels of randomly selected circRNAs were consistent with the preliminary RNA-seq data (Figure S2 and Figure 2C).

Figure 2 Differentially expressed circRNAs. (A) Hierarchical clustering of differentially expressed circRNAs between the AR group and control. The high and low relative expression levels of circRNAs were denoted by red and green respectively; (B) Volcano plot analysis based on circRNA expression levels. The red and blue points represent upregulated and downregulated circRNAs respectively (fold change ≥ 2.0, P value < 0.05). (C) The expression levels of the five upregulated and five downregulated circRNAs from the RNA-seq data were validated by qRT-PCR.

Figure S2 qRT-PCR validation of selected circRNAs. The relative expression of the ten circRNAs showed significance differences in AR compared to the control. ** indicates P < 0.01, *** indicates P < 0.001.

GO and KEGG pathways analyses

The circRNAs can exert their regulatory potency by interacting with their parental genes. Therefore, GO and KEGG analysis were performed to determine the biological functions of the parental genes. Further, 40 biological GO terms and 15 KEGG pathways were significantly enriched. Results of the GO analysis show that these parental genes were mainly enriched in terms associated with developmental process, cell junction, extracellular region, cytokine secretion, and signal transducer activity (Figure 3A). The most prominent GO terms in the biological process were the negative regulation of the T cell receptor signaling pathway (GO:0050860, P = 5.74E-4) and calcium ion transmembrane transport (GO:0070588, P = 1.62E-4), the most prominent GO terms in the cellular component were the endoplasmic reticulum membrane (GO:0005789, P = 1.60E-4), chromatin (GO:0000785, P = 4.76E-4), and membrane raft (GO:0045121, P = 1.97E-5) and the most prominent enriched GO terms in the molecular function were chromatin binding (GO:0003682, P = 4.27E-05) and calcium-release channel activity (GO:0015278, P = 5.01E-04). The biological pathways correlated with the parental genes were analyzed according to the KEGG database, and the top ranked five putative pathways were associated with the B cell receptor signaling pathway (nmu04662), calcium signaling pathway (nmu04020), MAP kinase (MAPK) signaling pathway (nmu04010), cyclic guanosine monophosphate/protein kinase G (cGMP- PKG) signaling pathway (nmu04022), and cholesterol metabolism (nmu04979) (Figure 3B). Some of these pathways have been reported to be involved in the pathogenesis of AR.

Figure 3 Enrichment analysis of circRNAs' parental gene function. (A) The biological process, molecular function, and cell components of these parental genes were annotated by using the GO database. (B) The biological pathways of these parental genes were annotated by using the KEGG database.

Construction of a circRNA–miRNA interaction network

Numerous experiments have reported that circRNAs can function as “miRNA sponges” to regulate the miRNA expression. The regulatory network of the top 10 upregulated and the top 10 downregulated circRNAs and their candidate microRNAs with good scores was explored in detail (Figure 4). Intriguingly, we found that some dysregulated circRNAs harbored several miRNA binding sites. For instance, the downregulated nmu_circ_0004141 (fold change = 4.14, P < 0.05) obtained higher binding scores with mmu-miR-135b-5p, mmu-miR-293-5p, mmu-miR-6769b-5p, mmu-miR-135a-5p, and mmu-miR-185-3p. Some of the predicted miRNAs, such as miR-135a, miR-143, and miR-466a, were conserved miRNAs both in murine and humans and have been previously demonstrated to function in AR development.^17–19 Therefore, the circRNAs related to these miRNAs may have an underlying role in AR development.

Figure 4 The view of a circRNA-miRNA interaction network. The network was comprised of 10 upregulated circRNAs (represented as red nodes) and 10 downregulated circRNAs (represented as blue nodes), the gray nodes represented the predicted miRNAs of the altered circRNAs.

Discussion

AR is characterized by the chronic inflammation of nasal mucosal caused by an IgE-mediated reaction. In this study, the mice in the AR group showed symptoms similar to the clinical behaviors, and their serum IgE levels were significantly higher than those of the healthy group, thereby justifying the appropriateness of our model. Previous studies revealed that the characteristic inflammation in AR is mainly induced by the predominance of the Th2 cell response and the inadequacy of the Th1 cell response.²⁰ Th2 cells can secrete cytokines, such as interleukin (IL)-4, IL-5, and IL-13. Furthermore, the Th2 cells induce an allergen-specific IgE production by triggering an IgE isotype shift in B cells, which subsequently lead to the maturation and aggregation of inflammatory cells following the release of histamine, leukotrienes, and prostaglandin.²¹ The suppression of Th2 has an enormous potential to be a candidate therapeutic strategy for treating AR.

To the best of our knowledge, this study is the first to perform RNA-seq analysis to identify the circRNA expression profiles of AR mice and controls. A total of 86 differentially expressed circRNAs were identified, of which 51 were upregulated and 35 were downregulated. The qRT-PCR detection of the randomly selected dysregulated circRNAs was consistent with the RNA-seq profiles, thereby highlighting the reliability of the sequencing data. Altered circRNAs were also found in all chromosomes, including the X and Y chromosomes. Therefore, all chromosomes were speculated to be involved in the abnormalities of AR, thereby suggesting that the pathogenesis of circRNAs in AR is complex and warrant further investigation.

GO and KEGG pathway analyses were conducted to explore the biological functions and vital pathways of circRNAs in AR. The GO enrichment showed that some parental genes were mainly involved in the transduction of calcium ion channels and cell signals.

Among the GO terms annotated in this study, T cell receptor signaling can act as a regulatory feedback machinery for T cell differentiation and homeostasis under differential cytokine circumstances.²² Among the KEGG pathways, B cell receptor signaling can function as a regulator of the immaturity and transition of B cells in vitro and play a role in the pathogenesis of autoimmunity.²³ Meanwhile, calcium signaling plays a paramount role in immunity. T cells generate a plethora of calcium-permeable channels at different locations with unique activation mechanisms that are essential for T cell activation, maturation, and production of cytokines.²⁴ Calcium signaling also plays a significant role in the development and fate of B cells.²⁵ Recent studies show that the cGMP-PKG pathway activated by antigen exposure plays a pivotal role in the augmentation of NO-induced vasodilatation of nasal mucosa, which in turn may contribute to the development of nasal obstruction in AR.²⁶ Taken together, these results indicate that circRNAs may play a regulatory function in AR.

Many studies have demonstrated the importance of circRNAs in the development and progression of diseases via their functions as miRNA sponges or competing endogenous RNA. Therefore, this research comprehensively profiled the interaction network of 20 candidates’ circRNAs and their potential binding miRNAs. Previous studies have also shown that multiple miRNAs were associated with the pathogenesis of AR, such as miR-16, miR-133b, and miR-375.^27–29 Some interacting microRNAs have been proven to play regulatory roles in balancing the Th1/Th2 immune condition. nmu_miR-135a is a predicated miRNA of mmu_circ_0004141 in this study, and previous studies reveal that miR-135a can efficiently correct Th1/Th2 imbalance by decreasing the expression of the GATA binding protein (GATA)-3 and IL-4 and by increasing the expression of T-bet and interferon- γ (IFN-γ) in AR mice.¹⁷ A similar study revealed that miR-466a-3p directly suppressed the transcription of GATA-3 and contributed to inadequate Th2 responses in nasal inflammation.¹⁹ miR-466a-3p is a target miRNA of mmu_circ_0009939, which is a downregulated circRNA identified in this study based on the miRNA binging site prediction. Therefore, it is evident that miRNAs play a critical role in the polarization and activation of T cells. The association of circRNAs and miRNAs indicate that the former may also have a regulatory role in AR. However, further research must be conducted to validate the circRNA–miRNA interaction.

Overall, this study is the first to screen the expression profiles of circRNA between AR murine models and the control. A series of aberrantly expressed circRNAs that may be closely associated with AR pathogenesis were identified. Further studies should be conducted to explore the underlying mechanisms of these circRNAs and to ascertain their potential as novel diagnostic biomarkers or therapeutic targets in AR.