The Transcriptome Characteristics of Severe Asthma From the Prospect of Co-Expressed Gene Modules

Rationale: Severe asthma is a heterogeneous disease with multiple molecular mechanisms. Gene expression studies of asthmatic bronchial epithelial cells have provided biological insights and underscored possible pathological mechanisms; however, the molecular basis in severe asthma is still poorly understood. Objective: The objective of this study was to identify the features of asthma and uncover the molecular basis of severe asthma in distinct molecular phenotype. Methods: The k-means clustering and differentially expressed genes (DEGs) were performed in 129 asthma individuals in the Severe Asthma Research Program. The DEG profiles were analyzed by weighted gene co-expression network analysis (WGCNA), and the expression value of each gene module in each individual was annotated by gene set variation analysis (GSVA). Results: Expression analysis defined five stable asthma subtype (AS): 1) Phagocytosis-Th2, 2) Normal-like, 3) Neutrophils, 4) Mucin-Th2, and 5) Interferon-Th1 and 15 co-expressed gene modules. “Phagocytosis-Th2” enriched for receptor-mediated endocytosis, upregulation of Toll-like receptor signal, and myeloid leukocyte activation. “Normal-like” is most similar to normal samples. “Mucin-Th2” preferentially expressed genes involved in O-glycan biosynthesis and unfolded protein response. “Interferon-Th1” displayed upregulation of genes that regulate networks involved in cell cycle, IFN gamma response, and CD8 TCR. The dysregulation of neural signal, REDOX, apoptosis, and O-glycan process were related to the severity of asthma. In non-TH2 subtype (Neutrophils and Interferon-Th1) with severe asthma individuals, the neural signals and IL26-related co-expression module were dysregulated more significantly compared to that in non-severe asthma. These data infer differences in the molecular evolution of asthma subtypes and identify opportunities for therapeutic development. Conclusions: Asthma is a heterogeneous disease. The co-expression analysis provides new insights into the biological mechanisms related to its phenotypes and the severity.


INTRODUCTION
Asthma is a chronic disorder, characterized by airway hyperresponsiveness (AHR) and remodeling with variable degrees of eosinophilic and neutrophilic inflammation resulting in significant morbidity and mortality (Wilson et al., 2006;Kim et al., 2010). It affects about 5% of the population (Global et al., 2017). According to the clinical characteristics, it is mainly divided into the acute and the non-acute asthma, which is further divided into mild, moderate, and severe asthma individual. About 5-10% of the patients do not respond well to standard treatment and have a poor prognosis (Higgins, 2003). The bronchial epithelial cells act as a physical barrier in airway immunity and as central modulators of inflammatory response (Hamilton et al., 2001). Environmental stimuli promote epithelial cell synthesis and secretion by a variety of mediators, such as cytokines, chemokines, reactive oxygen species, lipid, and peptide mediators and eventually involved in recruiting leukocytes, mucus secretion, vascular permeability, bronchoconstriction, and airway hyper-responsiveness (Whitsett, 2018).
Gene expression and genetic variation studies both indicate that asthma is a polygenic and heterogeneous disease with multiple molecular roots Belsky et al., 2013). Based on the gene profiles in bronchial epithelial cells associated with fractional exhaled nitric oxide (FeNO), Modena et al. identified five phenotypes of asthma. The results showed that a large number of individuals were severe asthma in each subtype (Modena et al., 2014). However, the typical characteristics of phenotype and the features related to severe asthma in phenotype were also unclear. Therefore, revealing these characteristics in each molecular subtype could be valuable for individualized treatment of severe asthma.
In recent years, the WGCNA ) (weighted gene co-expression network analysis) is a new system biology approach that can be used to identify co-expression gene sets that largely represent the typical biological characteristics in complex disease. Therefore, in this study, we used the coexpressed gene modules (GMs) that were used to uncover the typical features in subtype and severe asthma.

Study Population and Data Processing
As part of SARP (Severe Asthma Research Program), bronchial brushing samples and matching demographic data were obtained from 155 participants (129 asthmatics and 26 healthy subjects) from 2009 to 2011. Gene expression of the SARP and external cohorts are available online (GEO database; http://www.ncbi.\.

K-Means and Limma Analysis
According to the transcription of bronchial epithelial cells, the k-means method integrated in the ConsensusClusterPlus (Wilkerson and Hayes, 2010) package was adopted to identify stable subtypes, and the stability was evaluated by iterating for 1,000 times at a sub-sampling rate of 0.95. A total of 4,650 (MAD: median absolute deviation >0.5) genes were used as input. Starting from k 4-5, a significant improvement in clustering stability can be observed, but it has no effect on k > 5 ( Figure 1C). They were termed by asthma subtype. After the establishment of these ASs, the Bayesian method in Limma package was used to select the differentially expressed genes (DEGs) between each ASs and the normal. The cutoff setting: FC ≥ log (1.5), FDR ≤ 0.05.

WGCNA Co-expressed Analysis
Using the default parameter setting and the 2,664 DEGs selected in ASs, the WGCNA was performed. This method clusters genes into modules using a topological overlap measure (TOM) . The TOM was a highly robust interconnection measurement method that essentially provided a measure of the connection strength between two adjacent genes and all other genes in a network. Genes were clustered using 1-TOM as the distance measure and GMs were defined as branches of the resulting cluster tree using a dynamic branch-cutting algorithm. Based on the dysregulated direction of each gene in asthma, the genes in each co-expressed module were split into 2 GMs.

Gene Set Variation Analysis
Gene Set Variation Analysis (GSVA) was performed using the R package "GSVA" (Hänzelmann et al., 2013) (function gsvaarguments: method "gsva", mx. diff TRUE). GSVA implements a non-parametric unsupervised method of gene set enrichment that allowed an assessment of the relative enrichment of a selected pathway across the sample space. The output of GSVA was a gene set by sample matrix of GSVA enrichment scores that were approximately normally distributed.
GSVA enrichment scores were generated for each gene set using the normalized gene expression data.

Pathway Enrichment Analysis
Using the clusterProfiler (Yu et al., 2012) package, the functional enrichment analysis was performed for the up-and downregulated GMs. Significance cutoff was defined as FDR <0.05 for multiple testing.

Correlation Analysis
The correlation analysis was performed by Pearson. The GSVA score in each gene set represented its overall expression in individual. The C2 (curated gene sets) dataset was downloaded from Molecular Signatures Database (MSigDB), which is a collection of annotated gene sets for pathway analysis (http:// software.broadinstitute.org/gsea/msigdb).

The Co-expression Features of AS
A total of 2,664 DEGs were detected among these five ASs compared to the normal samples. Using 2,664 DEG as input, a FIGURE 1 | The heatmap of molecular subtype in asthma and the percentage of severe asthma in each cluster. The heatmap of molecular subtype in asthma and the percentage of severe asthma in each cluster (A) The k-means identified five stable asthma subtypes (AS1-5), the heatmap was derived by 2,664 DEGs selected in each subtype. To define the Th2 subtype, the unsupervised hierarchical clustering was performed based on the microarray expression levels of periostin (POSTN), channel regulator 1 (CLCA1), and serpin peptidase inhibitor, clade B, member 2 (SERPINB2). We named these subtypes: (1) Th-H, High; (2) Th-M, Moderate; (3) Th-L, Low. According to the severity of the disease, the samples were divided into five groups: Normal for normal group; Mid/Mod-noICS for mild and moderate without ICS treatment; Mid + ICS for mild and ICS treatment group; Mod + ICS for moderate and ICS treatment group; and Severe for severe asthma group (B) The percentage of severe asthma in each cluster (red, normal samples, green: non-severe asthma individuals; blue: severe asthma individuals).
Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 765400 total of 15 coordinately expressed GMs representing distinct biological processes were obtained. In these 15 GMs with only upregulated genes, 10 discriminated these five asthma clusters ( Figure 2). In the validation dataset, compared to the normal group, the overall expression level of these modules was consistent with the trend in this study.

Normal-like Subtype: AS2
The Normal-like subtype was the most similar to the normal group, with the mildest clinical symptoms and the fewest differential expressed genes. The top three upregulated genes were PHACTR3 (Itoh et al., 2014), SLCO1B3, and GNMT, which were related to the response of glucocorticoids.

Neutrophils-type Subtype: AS3
The GM10 was a typical feature of Neutrophils-Type. This module only contained eight upregulated genes, including SLCO1B3, PHACTR3, TPO, and FKBP5 (Binder, 2009), which were also related to the glucocorticoids and severity. The TPO was a marker that related to the severity of asthma (Voraphani et al., 2014).

The Characteristics Related to Severity in Asthma
In GMs with upregulated genes, 9 GMs (0, 1, 3, 4, 6, 8, 10, 12, and 14) were positively correlated to the severity and positive association with the use of ICS and OCS. They were also the FIGURE 2 | The heatmap of modules with upregulated genes in ASs. The overall expression is represented by red and blue; the red indicates high expressed in cluster, and the blue indicates low expressed in cluster. M0-M14 represents 15 co-expressed modules with upregulated genes.
The comparison between mild-mod-ICS and mild-mod-noICS showed that ICS/OCS significantly reduced the expression of 3 GMs (2, 7, and 13) with upregulated genes, while increasing the expression of GM4 with upregulated modules.
In "Neutrophils-Type" and "Interferon-Th1" (non-Th2) phenotypes, two specific modules (GM13-up and GM12-down) show different expression between severe asthma and non-severe asthma. The GM13-up related to interleukin-26 (IL26, PIK3R5, and LRRC2) was much higher expressed in severe individuals while the GM12-down module related to the nervous system (10%, p 1.77E-FIGURE 5 | Seven modules expression related to asthma severity. The geometric means were measured according to asthma severity. These classes included healthy controls (Normal), mild-to-moderate asthma not on inhaled corticosteroids (Mild-Mod-noICS), mild-to-moderate asthma on inhaled corticosteroids (Mild-Mod-ICS), and severe asthma (SA).

DISCUSSION
Asthma is a heterogeneous disease with multiple immune and non-immune mechanisms (Ito et al., 2004;McKinley et al., 2008). This study shows that the transcriptome of bronchial epithelial cells was related to ASs and the severity.
Using cluster analysis, five stable ASs were obtained. Multiple clinical features had significant differences among these ASs. These results were similar to that of previous studies . Using the expression of 2,664 DEG profiles as input, the WGCNA co-expression analysis was performed and the module was split based on the dysregulated direction of each gene; 30 GMs were obtained, 15 of which contain only upregulated genes and the other 15 contain downregulated genes ( Table 3). These split modules were essential for the description of the typical characteristics in ASs. In each co-expressed module, there was a negative correlation between the two gene sets with up-and FIGURE 6 | Six modules expression in each phenotype between severe and non-severe groups. The geometric means were measured according to each cluster. These classes included healthy controls (Normal), five asthma phenotypes (AS1-5), non-severe asthma, and severe asthma (SA).
Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 765400 8 downregulated genes, respectively. For example, in GM9 with upregulated genes, multiple Th2-related marker genes (CLCA1, POSTN, etc.) were included, while in GM9 with downregulated genes, the MUC5B (Zhang et al., 2019), SLC28A3, and CSGALNACT1 genes were included and closely related to the reduction of airway defense response (Ridley and Thornton, 2018;Rojas et al., 2019). The MUC5B plays a key role in airway defense. The lack of MUC5B leads to lung inflammation, impaired immune balance, and chronic infection mediated by a variety of bacteria (Roy et al., 2014). These two aspects (up and down features) might be two effective strategies for individualized treatment in asthma.
Among these 15 modules containing upregulated genes, 4 GMs (2, 3, 5, and 7) were typical immune-related modules and had obvious subtype distribution characteristics. The GM2 and GM3 were typical characteristics of "Phagocytosis-Th2" while the GM5 and GM7 were the typical features of "Interferon-Th1". The function enrichment analysis showed that the GM2 and GM3 were mainly related to the receptor-mediated endocytosis and antigen presentation via MHC-II. The GM5 and GM7 were highly expressed in "Interferon-Th1" and mainly related to type I interferon response and T-cell toxicity. Although the "Mucin-Th2" was a typical Th2 phenotype, the four immune-related modules mentioned above were not significantly increased in this cluster, while the increased glycosylation of O-type glycans, amino sugar and nucleoside sugar metabolism, proteolysis, and unfolded protein reaction were highly expressed in this cluster and positively correlated to the expression of Th2 markers. It suggested that these biological functions could be coordinated with the Th2 signals and related to the physio-pathological mechanisms in "Mucin-Th2". The increased mucus in airway was a typical clinical feature of Th2 asthma (Dunican et al., 2018;Lambrecht et al., 2019). Studies had shown that the galectin-10 (Galectin-10 and Gal10) crystal structure in airway mucus stimulates the immune system and induces the changes in airway inflammation and mucus secretion (Nyenhuis et al., 2019;Persson et al., 2019). Clearing the crystal structure can effectively relieve airway inflammation and asthma symptoms.
Correlation analysis showed that the apoptotic process and O-glycan processing were positively correlated to asthma severity, while the cell adhesion, innate immune response, potassium and sodium ion transmembrane transport, cellular oxidant detoxification, and neuron signal were negatively correlated to asthma severity. Correcting the imbalance of oxidation and antioxidation in the lung may be an important method to relieve asthma symptoms in clinic. GSH is the most important antioxidant in lung tissues (Brigelius-Flohé and Maiorino, 2013). GSH can inhibit a variety of pathogen replication and survival, and increasing the GSH can effectively improve the body's ability to resist foreign microorganisms (Fitzpatrick et al., 2012). The inhibition of the activity of the CYP450 pathway could destroy the phagocytosis of macrophage and reduces the clearance efficiency of inflammatory stimuli (Bystrom et al., 2013).  In "Phagocytosis-Th2", the GM10-down was differentially expressed between severe and non-severe groups. In "Mucin-Th2", the O-glycan processing (GM1-Up), apoptotic process (GM8-Up), and oxidation-reduction process (GM8-down) were significantly different between severe and non-severe groups. The oxygen free radical increase in bronchoalveolar lavage fluid and peripheral blood associated with the severity of disease (Mossberg et al., 2009;Sangiuolo et al., 2015), especially in a typical Th2 phenotype (Huang et al., 2019).
In "Neutrophils-Type" and "Interferon-Th1" (non-Th2) phenotype, interleukin-26 (IL26)-related function was upregulated and related to the severity of asthma. IL-26 is a member of IL-10 cytokine family, is abundant in human airways, and induces the production of pro-inflammatory cytokines (Louhaichi et al., 2020). Stimulation of cultured CD4 + T cells with monocyte by recombining IL-26 promoted the generation of RORγ Th17 + cells, inducing the production of IL-17A, IL-1β, IL-6, and TNF-α (Louhaichi et al., 2020). Therefore, IL-26 could appear as a novel pro-inflammatory cytokine, produced in airways, and may be a promising target to treat inflammatory asthma.
Although we clustered the transcriptome of bronchial epithelial cells and revealed the typical features in five stable subtypes, the heterogeneity in asthma is much higher than that in subtypes. The characteristics in individuals were more likely a mixture of typical features in multiple subtypes. For asthma, distinguishing subtypes was only a powerful method for uncovering the heterogeneity of complex diseases. Therefore, the individualized analysis based on phenotypes in asthma was a powerful tool for individualized diagnosis and treatment.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article.

AUTHOR CONTRIBUTIONS
Acquisition of data: BL accessed to the transcriptom data online; Conception and design: YZ, X-FJ, BL; Analysis and interpretation: BL, W-YZ, W-XS, YZ, LQ, Y-MH, W-QL, DL; Wrote the article: BL, W-XS; Approved and edited the article: All authors approved the article.

FUNDING
This study is funded by the National Natureal science foundation number 81171657, No,30371364.