Bioinformatics analysis of mRNA profiles and identification of microRNA-mRNA network in CD4+ T cells in seasonal allergic rhinitis.

OBJECTIVE: We aimed to discover potential circulating genes and non-coding molecules (micro RNA [miRNA] and circular RNA [circRNA]) in CD4+ T cells in relation to seasonal allergic rhinitis (SAR).
METHODS: Microarray data of GSE50223 were obtained from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) during and outside the pollen season were analyzed using R software and by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathway analyses. The protein-protein interactions, modules, miRNAs targeting DEGs, merged miRNA-DEG networks, and circRNAs targeted with miRNAs were further analyzed.
RESULTS: We identified 211 DEGs during the pollen season and eight DEGs outside the season, of which only MMP12, NR4A2, and CD69 were differentially expressed both during and outside the pollen season. DEGs during the pollen season were enriched in the GO categories 'neutrophil degranulation', 'neutrophil activation involved in immune response', 'neutrophil mediated immunity', and 'neutrophil activation'. A significant module was identified with key nodes of CDK6 and hsa-miR-29b-3p. Six significant circRNAs were also identified.
CONCLUSIONS: Some genes, miRNAs, and circRNAs in CD4+ T may play vital roles in SAR and may thus be potential targets for the prevention and treatment of SAR.

Introduction

Allergic rhinitis (AR) is an IgE‐mediated type-2 upper airway inflammatory disorder with increasing prevalence worldwide, which is characterized by nasal congestion, itching, sneezing, and rhinorhea.[1-3] AR can be divided into seasonal AR (SAR) and perennial AR according to allergen tests. The symptoms are commonly controlled by environmental control measures, intranasal corticosteroids, and antihistamines, while allergen immunotherapy is currently the only disease‐modifying treatment option for AR patients.[4-6] However, AR confers a significant healthcare burden. There is thus a need to investigate the mechanisms and identify better biomarkers for AR to broaden the available therapeutic approaches.

CD4+ T cells play a central role in atopic diseases including AR, allergic asthma, and allergic dermatitis.[8-10] CD4+ T cells can be subdivided into effector cells and regulatory T (Treg) cells, based on their function. These subsets perform different functions, with effector cells mainly providing protection against exogenous offending agents, and Treg cells being involved in avoiding autoimmune reactions and stopping the effector response against exogenous antigens. Human effector CD4+ T lymphocytes can also be classified into lineages according to their immunological functions, based on their cytokine, transcription factor, and homing receptor expression profiles, including T helper (Th) lymphocytes, Th1 and Th2 cells, Th17 subset, Th22, Th9, and T follicular helper lymphocytes. In addition to their protective functions, CD4+ T cells also mediate and participate in autoimmune and allergic diseases and a variety of inflammatory diseases, and allergen-specific Th2 lymphocytes play a crucial role in allergic inflammation.[8,9]

Th2 cytokines, including IL-4, IL-5, and IL-13, have been found to be essential factors in the pathogenesis of AR.[1,3] Although biologicals, such as drugs targeting IL-5, have been confirmed to be a useful therapeutic strategy for severe asthma, their use in AR is restricted.[11,12] However, IL-35 was proven to inhibit allergic responses and its intranasal administration improved AR symptoms in murine models.[13,14]

MicroRNAs (miRNAs) are small non-protein-coding RNA molecules that can decrease gene expression by degrading mRNA or inhibiting translation through binding to the 3′-untranslated region of their target mRNAs. Has-let-7 family members were shown to be decreased in both animal models and human studies of asthma and AR. miR-487b inhibited the IL-33/ST2 signaling pathway and miR-133b inhibited Nlrp3 inflammasome-mediated inflammation, and both were reported to mitigate AR in an ovalbumin-induced animal model.

Although previous studies have reported some pathomechanisms of AR, several aspects remain unclear. RNA sequencing is a powerful and unbiased high-throughput method for exploring candidate transcripts. In this study, we downloaded the GSE50223 gene expression profile dataset for CD4+ T cells from SAR samples and healthy controls from the Gene Expression Omnibus (GEO) database and analyzed the data using bioinformatics analysis. Key genes and their targeted non-coding RNAs were identified as possible biomarkers of SAR, thus improving our understanding of the genetic mechanisms of the disease.

Methods

Data collection and screening of differentially expressed genes (DEGs)

The gene expression profiling dataset (GSE50223) was obtained from the National Center for Biotechnology Information GEO database (https://www.ncbi.nlm.nih.gov/geo), based on the GPL6884 platform (Human WG-6 v3.0 expression Beadchip; Illumina, San Diego, CA, USA). Using this dataset, we obtained data for peripheral blood mononuclear cells from 21 healthy controls (H) and 21 untreated SAR patients (P) outside the pollen season, challenged using diluent (phosphate-buffered saline) or allergen (grass pollen extract) for 7 days. CD4+ T cells were then purified for gene expression profiling. DEGs between patients with challenge allergen (PA) and healthy controls with challenge allergen (HA), and between patients with challenge diluent (PD) and healthy controls with challenge diluent (HD), were identified using the limma package in R software. Significant differences were determined by a P-value of <0.05 and fold-change (FC) ≥1.5 or ≤2/3. Subjects with challenge diluent (PD and HD) reflected the status of subjects out of the season, and subjects with challenge allergen (PA and HA) reflected the status during the pollen season.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway assessments

We further explored the functions of the DEGs by GO and KEGG pathway analyses using the cluster profile package in R software. GO enrichment analysis included biological process (BP), molecular function (MF), and cellular component (CC). The KEGG pathway database records networks of cellular molecular interactions. The downloaded results of the GO and KEGG analyses were based on an adjusted P value <0.05, which was considered to be significant.

Protein–protein interaction (PPI) network construction and module selection

DEGs between PD and HD were rare and we therefore only analyzed DEGs between PA and HA to show the PPI relationship. Corresponding DEGs were uploaded to the Search Tool for Retrieval of Interacting Genes (STRING) database (version 11.0; https://string-db.org) to construct functional protein association networks. The validated connection was confirmed by an interaction score of genes >0.4. PPI networks were visualized in Cytoscape (v.3.7.1) and DEG modules were distinguished using the Molecular Complex Detection (MCODE) plug‑in.[21,22] The criteria were set as scores >4 and nodes >4.

Prediction of potential DEG target miRNAs and construction of miRNA–mRNA regulatory network

Based on DEGs between HA with PA, we predicted the target miRNAs of the upregulated and downregulated genes using Targetscan Human 7.1 (http://www.targetscan.org/), and constructed a miRNA–mRNA interaction graph using Cytoscape (v.3.7.1) to show the interactions intuitively. PPI pairs and miRNA–DEG pairs were then merged to create a miRNA–DEG network, which was also visualized using Cytoscape software (v.3.7.1). Module clustering analysis for the network was then performed using the MCODE plug-in[21,22] to detect potential functional modules in the network.

Prediction of potential circular RNAs (circRNAs)

The CircBank database (http://www.circbank.cn) is a comprehensive database of human circRNAs. We screened CircBank using previous screening criteria to determine the circRNAs responsible for the regulation of the target miRNAs.

Results

DEG screening between SAR patients and healthy controls

Only eight DEGs, including three upregulated and five downregulated genes, were selected between the PD and HD groups (Figure 1a). However, 211 DEGs between the PA and HA groups were identified from the GSE50223 dataset, including 163 upregulated and 48 downregulated DEGs (Figure 1b). Among these, CD14 (encoding CD14) was the most significantly upregulated gene (with the highest logFC, 2.14; P  =  1.81E-06) and IL17RB (interleukin [IL]-17 receptor B) was the most significantly downregulated gene (with the smallest logFC, −2.43; P = 3.00E-13). We also compared the DEGs in each pairwise analysis. NR4A2 (nuclear receptor subfamily 4 group A member 2, NR4A2) and CD69 (CD69) were downregulated in both HD vs PD and HA vs PA subjects, whereas only MMP12 (encoding matrix metalloproteinase [MMP] 12) was upregulated (Figure 1c, 1d).

Figure 1.

Alterations in mRNA expression profiles in patients with seasonal allergic rhinitis. Volcano plots of differentially expressed genes (DEGs) between (a) patients with challenge diluent (PD) and healthy controls with challenge diluent (HD), and (b) between patients with challenge allergen (PA) and healthy controls with challenge allergen (HA). (c) Upregulated and (d) downregulated DEGs in HD-PD vs HA-PA. Red dots indicate upregulated DEGs and blue dots indicate downregulated DEGs.

Functional and pathway enrichment analyses

After uploading the DEGs to the cluster profile package in R software, we identified the GO categories with P < 0.05. The significant BP terms for HD vs PD and the top 20 significant BP terms for HA vs PA were selected separately from the corresponding GO categories (Table 1). GO analysis revealed that HD vs PD DEGs were significantly associated with ‘defense response’, ‘response to external stimulus’, ‘immune response’, ‘response to stress’, and ‘regulation of multicellular organismal process’ and HA vs PA DEGs were significantly enriched in ‘immune response’, ‘immune system process’, ‘myeloid leukocyte activation’, ‘leukocyte activation’, ‘neutrophil degranulation’, ‘neutrophil activation involved in immune response’, ‘neutrophil-mediated immunity’, ‘neutrophil activation’, ‘leukocyte degranulation’, ‘granulocyte activation’, ‘cell activation’, ‘myeloid cell activation involved in immune response’, ‘myeloid leukocyte-mediated immunity’, ’immune effector process’, ‘secretion’, ‘secretion by cell’, ‘leukocyte activation involved in immune response’, ‘cell activation involved in immune response’, ‘leukocyte mediated immunity’, and ’defense response’.

Table 1.

Top enriched Gene Ontology (GO) terms among differentially expressed genes in patients with challenge diluent (PD) vs healthy controls with challenge diluent (HD) and patients with challenge allergen (PA) vs healthy controls with challenge allergen (HA).

GO ID	Description	Gene count	P-value
Top 5 enriched GO terms among PD-HD DEGs
GO:0006952	defense response	4	0.008
GO:0009605	response to external stimulus	4	0.017
GO:0006955	immune response	4	0.021
GO:0006950	response to stress	5	0.022
GO:0051239	regulation of multicellular organismal process	4	0.038
Top 20 enriched GO terms among PA-HA DEGs
GO:0006955	immune response	107	<4.44E-35
GO:0002376	immune system process	121	2.65E-34
GO:0002274	myeloid leukocyte activation	55	1.84E-30
GO:0045321	leukocyte activation	72	1.35E-28
GO:0043312	neutrophil degranulation	47	1.55E-28
GO:0002283	neutrophil activation involved in immune response	47	2.45E-28
GO:0002446	neutrophil mediated immunity	47	4.48E-28
GO:0042119	neutrophil activation	47	6.03E-28
GO:0043299	leukocyte degranulation	48	8.71E-28
GO:0036230	granulocyte activation	47	1.17E-27
GO:0001775	cell activation	74	3.20E-27
GO:0002275	myeloid cell activation involved in immune response	47	1.47E-26
GO:0002444	myeloid leukocyte mediated immunity	47	2.35E-26
GO:0002252	immune effector process	69	6.49E-26
GO:0046903	secretion	75	8.66E-26
GO:0032940	secretion by cell	72	1.04E-25
GO:0002366	leukocyte activation involved in immune response	51	1.33E-25
GO:0002263	cell activation involved in immune response	51	1.66E-25
GO:0002443	leukocyte mediated immunity	58	1.13E-24
GO:0006952	defense response	81	1.44E-24

In addition, the enriched pathways for HD vs PD and the top 20 significant enriched pathways for HA vs PA were selected separately from the significant KEGG pathways (Figure 2a, 2b). HD vs PD DEGs were significantly enriched in ‘Aldosterone synthesis and secretion’, ‘Asthma’, ‘Parathyroid hormone synthesis, secretion and action’, ’Toll-like receptor signaling pathway’, and ‘Lysosome’, and HA vs PA DEGs were associated with ‘Cytokine-cytokine receptor interaction’, ‘Osteoclast differentiation’, ‘Complement and coagulation cascades’, ‘Lysosome’, ‘Transcriptional misregulation in cancers’, ‘Tuberculosis’, ‘IL-17 signaling pathway’, ‘Chagas disease’, ‘Acute myeloid leukemia’, ‘Measles’, ‘Salmonella infection’, ‘Cholesterol metabolism’, ‘Malaria’, ‘Rheumatoid arthritis’, ‘Toll-like receptor signaling pathway’, ‘Leishmaniasis’, ‘Pertussis’, ‘Staphylococcus aureus infection’, ‘African trypanosomiasis’, and ‘HIF-1 signaling pathway’. The KEGG pathway annotations are indicated in Figure 2c and 2d.

Figure 2.

Significant pathways identified in the Kyoto Encyclopedia of Genes and genomes (KEGG) database. (a) The most significant enriched pathways in healthy controls with challenge diluent (HD) vs patients with challenge diluent (PD). (b) KEGG pathway annotation of (a). (c) Top 20 enriched pathways in healthy controls with challenge allergen (HA) vs patients with challenge allergen (PA) and (d) KEGG pathway annotation of (c).

PPI network construction and module selection

A PPI network of the HA vs PA DEGs was produced using the STRING online tool with a cutoff score >0.4. A total of 1300 pairs of PPIs were identified and a big pairing picture was generated using Cytoscape (Figure 3a). MCODE analysis revealed a total of 10 available modules, of which the top five significant modules were selected (Figure 3b–f).

Figure 3.

Protein–protein interaction regulatory network of differentially expressed genes and top five modules in seasonal allergic rhinitis. (a) Protein–protein interaction regulatory network. (b) Module 1, (c) Module 2, (d) Module 3, (e) Module 4, and (f) Module 5. Nodes with higher degree values depicted by larger circles. Edges/lines indicate regulatory association between any two nodes. Red and green nodes represent upregulated and downregulated genes, respectively.

Construction of miRNA–DEG regulatory network

Among the 163 upregulated and 48 downregulated HA vs PA DEGs, 313 and 218 respective miRNAs were predicted by Targetscan, 1493 and 554 of the identified miRNA–DEG pairs were visualized using Cytoscape, and two big pictures were generated (Figure 4).

Figure 4.

Regulatory network of differentially expressed genes (DEGs) and their binding microRNAs (miRNAs) in seasonal allergic rhinitis. (a) Upregulated and (b) downregulated DEGs and their binding miRNAs. Red and green nodes represent upregulated and downregulated DEGs, respectively; blue and yellow nodes represent binding miRNA of upregulated and downregulated genes, respectively.

Key functional-module accession and potential circRNA predictions

Given that the above networks contained extremely large amounts of complex information, we performed module clustering based on the miRNA–DEG and PPI networks using the MCODE plugin in Cytoscape to assess the key functional modules. One module was identified and showed upregulated DEGs including MMP9, CD68, S100A8 (encoding S100 calcium-binding protein A8), IL1B (IL-1β), and TLR8 (Toll-like receptor 8), downregulated DEGs including CDK6 (cyclin-dependent kinase 6, CDK6), IFNG (interferon-γ, IFNG), and EOMES (eomesodermin, EOMES), and targeted miRNAs including hsa-miR-29b-3p and hsa-let-7 family (including hsa-let-7a-5p, hsa-let-7b-5p, hsa-let-7c-5p, hsa-let-7d-5p, hsa-let-7e-5p, hsa-let-7f-5p, hsa-let-7g-5p, and hsa-let-7i-5p) (Figure 5a). The distributions of the above eight DEGs among the top 20 significant BP terms for HA vs PA GO categories are shown in Figure 5b. Furthermore, we determined the circRNAs responsible for the regulation of hsa-miR-29b-3p by screening CircBank using a detailed screening system, and finally screened out the six most promising candidates, including hsa_circ_0095386, hsa_circ_0044548, hsa_circ_0081146, hsa_circ_0057377, hsa_circ_0057450, and hsa_circ_0056444 (Figure 5c).

Figure 5.

Key functional module accession and potential circRNA prediction. (a) Most significant module from the intact microRNA–differentially expressed gene (DEG) network. (b) Distribution of the eight DEGs in the top 20 significant Gene Ontology (GO) Biological Process terms and (c) Schematic of screening rules for circular RNAs (circRNAs) that might regulate miR-29b-3p.

Discussion

AR is a common allergic upper respiratory disease involving multiple pathogenic elements.[1,3] Patients with AR share common pathophysiologies, including heightened nasal mucosa hyperresponsiveness and heightened reactivity to various stimuli.[1,2] However, most previous research has focused on analyzing individual or limited numbers of markers,[14,18,26] and whole-genome approaches using well-powered sample sizes are necessary to elucidate the gene signatures of AR. The findings of the current study may lead to an improved and more comprehensive understanding of the genetic mechanisms of circulating CD4+ T cells in relation to AR.

The present study showed that SAR was associated with different gene expression profiles during and outside the pollen season, and only MMP12, NR4A2, and CD69 were expressed differentially between SAR patients and controls both during and outside the pollen season. MMP12, a member of the MMP family, is structurally related to extracellular matrix-degrading enzymes and is thought to be involved in the tissue remodeling process and the development of inflammatory responses. Although it has not been investigated in relation to AR, MMP12 was suggested as a potential therapeutic target for the treatment of allergic contact dermatitis. NR4A2 is part of the nuclear orphan receptor 4a family. Knockdown of NR4A2 suppressed the IgE-induced upregulation of tumor necrosis factor-α mRNA expression in mast cells, and downregulation of NR4A2 in CD4+ T cells has been demonstrated in atopic dermatitis. CD69 is involved in lymphocyte proliferation and acts as a signal-transmitting receptor in lymphocytes, natural killer cells, and platelets, and has been reported to limit allergic asthma and skin contact hypersensitivity. MMP12, NR4A2, and CD69 might thus play important roles in the development of SAR. We accordingly showed separation of patients and controls by these three genes both during and outside the pollen season. The “allergic march” has become a hot research topic in recent years, leading to improved overall understanding of this allergic disease. However, the mechanism of the allergic march is not completely clear. Interestingly, the three DEGs (MMP12, CD69, and NR4A2) that were identified both in and out of the pollen season in patients with SAR are also associated with dermatitis,[31,33,34] thus supporting their participation in the mechanism of the allergic march. However, further in-depth investigations are needed to clarify the mechanism underlying the correlations between MMP12, CD69, NR4A2 and the allergic march.

Patients with AR demonstrate comparable systemic and local eosinophil inflammation, but relatively little is known about the contribution of neutrophil inflammation in AR. Neutrophil–lymphocyte ratios were recently found to be higher in patients with moderate/severe AR, and might thus be a useful marker for the diagnosis of persistent AR.[36,37] In the current study, the results of GO and KEGG analyses reflected the functional and pathway information for 211 DEGs during the pollen season, with ‘neutrophil degranulation’, ‘neutrophil activation involved in immune response’, ‘neutrophil mediated immunity’, and ‘neutrophil activation’ among the top 20 GO analysis terms, indicating the importance of CD4+ T cell-mediated neutrophil biological processes in the development of SAR, and providing novel ideas for understanding the mechanisms of this disease. Neutrophil inflammation might play an important role in the potential induction of immune intolerance associated with persistent airway inflammation in SAR.

The above combined network analysis showed that all eight genes in the significant module were enriched in the top 20 GO analysis. MCODE analysis also showed that CDK6 and hsa-miR-29b-3p might be the core of the module, and hsa-miR-29b-3p was associated with the CDK6, IFNG, and EOMES genes. CDK6 is a classic cell cycle kinase that forms complexes with D-type cyclins to regulate cell transition from G1 to S phase. CDK6 has been shown to regulate several genes, including those encoding early growth response-1 and the angiogenic factor vascular endothelial growth factor A. IFNG bridges the innate and adaptive immune system and is a key cytokine in the development of a Th1-type response, and also plays a role in the regulation of local leukocyte–endothelial interactions. Serum IFNG was negatively correlated with nasal symptom scores. Together with the genes for IFNG, T-box expressed in T cells, and IL-12 receptor β1, EOMES is a prototypic Th1 gene, and suppression of EOMES can reduce the development of Th2 cells and allergic airway inflammation. In IgA nephropathy, downregulation of miR-29b-3p caused by CDK6 overexpression was proven to promote nuclear factor-κB signaling by phosphorylating p65 and enhancing inflammation, while downregulation of miR-29b-3p also activated CDK6 and inhibited cell-cycle arrest at G1 phase in breast cancer cells. All these results suggest that the key nodes of CDK6 and hsa-miR-29b-3p might play vital roles as mediators of inflammation reactions in SAR.

Recent evidence has shown that circRNAs can act as miRNA sponges, binding to miRNAs to inhibit their functions in many types of human diseases, such as diabetes, cancer, and osteoarthritis. Overexpression of circDdx17 inhibited the expression of miR-17-5p and alleviated the allergic condition in an AR animal model. However, the contribution of circRNAs in AR remains largely unknown. The six most promising candidate circRNAs screened out in the current study might thus lead to further investigations of their roles in the pathogenesis of SAR.

The present study explored the mechanisms underlying seasonal AR through various bioinformatics approaches, to lay theoretical foundations for further research into promising targets for SAR therapy. However, this study had some limitations, including the lack of experimental trials. Both in vivo and in vitro follow-up studies are therefore needed to support our results, and additional investigations into the molecular mechanisms of SAR pathogenesis are needed to assess the functional value of these results.

Conclusion

A total of 211 circulating DEGs were identified between patients with SAR and controls during the pollen season, and only MMP12, NR4A2, and CD69 were differentially expressed both during and outside the pollen season. A significant module was identified with key nodes of CDK6 and hsa-miR-29b-3p, suggesting that these might play vital roles in mediating inflammation reactions in SAR. Moreover, hsa-miR-29b-3p might also be implicated in SAR through regulating CDK6, IFNG, and EOMES. In addition, the circRNAs identified in this study may be candidate genes for subsequent experimental research.