Integration of Airway Inflammation and Remodeling Mechanisms Specific to Eosinophilic Asthma Through Differential Co-Expression of Genes in Bronchial Brush Biopsy Samples


 Background

Heterogeneity of asthma complicates search for targeted treatment against airway hyperresponsiveness and remodeling. We conducted a systems biology approach study to establish differential co-expression of genes (DCG) in eosinophilic and non-eosinophilic asthma patients and infer their role in the disease.
Materials and Methods

N = 40 Caucasian adult moderate to severe non-smoking asthma patients (half with eosinophilic asthma) undergone bronchial brush biopsy sampling for mRNA expression using hybridization to cDNA microarray. Mechanistic interpretation of DCG was inferred from existing literature.
Results

Differentially co-expressed genes bear significance in airway viral infection (ATP1B1, EPS15), arachidonic acid metabolism (CLC, FADS6), cell migration (EPS8L1, STOML3, RhoBTB2), surface receptors endocytosis (STRN4, EPS15, ATP1B1) or decreased expression (CCT7), oxidative stress (DIO3, RhoBTB2), decreased adhesion (ATP1B1, RAPH1, STOML3), epithelial-mesenchymal transition (ASB3, RADX, CCT7, MRPL14, PPP2R3B, RPS13, SLC19A1), myofibroblast differentiation (CCT7), smooth muscle proliferation (ASB3, ATP1B1), airway hyperreactivity (RECK, STOML3, ATP1B1, OR52I1), extracellular matrix remodeling (FBN3, RECK), angiogenesis (GPI, RhoBTB2) and neuronal pathogenesis of asthma (OR52I1, STRN4, TTC3P1, GPI, CABP5) and were linked to asthma in genome- (MRPL14, ASB3, RPS13) and epigenome-wide (CLC, EPS15, GPI, SSCRB4, STRN4) association studies. Signaling pathways involved (especially TGF-β/Smad2/3) are inferred from the co-expression pattern.
Conclusion

Activity of genes and pathways of known or tentative role in asthma pathogenesis was established in regard to a condition cognizable in clinical practice.

subepithelial brosis, epithelial-mesenchymal and broblast-to-myo broblast transition, hyperplasia of smooth muscle and goblet cells, neovascularization and cartilage degeneration 1 , leading to hyperresponsiveness and persistent air ow limitation. Current treatment options have limited in uence on remodeling and do not target it directly as mechanisms underlying the process are still investigated 23 . Since no accurate biomarker or clinical predictor of ongoing remodeling was agreed upon to date, there are no methods of selecting patients who could bene t the most from future experimental e.g. biological therapies before the process becomes clinically symptomatic or apparent in imaging studies 4 .
As eosinophilic asthma (EA) accounts for approximately half of the severe disease cases 5 , search for therapeutic targets behind airway in ammation and remodeling speci c to this phenotype presents an important scienti c challenge.
In the present hypothesis-generating study we utilized systems biology approach to identify genes with coordinated expression in the airways of EA and non-EA patients, established networks of their binary interactions and regulation using curated resources and conducted a systematic review of the relevant literature to infer their possible role in the disease pathogenesis in search for novel regulators, biomarkers and mechanisms of airway in ammation and remodeling.

.1 Patient selection
The study included 40 non-smoking patients aged 20-70 years with at least 10-year history of asthma diagnosed by a physician in accordance with Global Initiative for Asthma (GINA) 2018 guidelines. The exclusion criteria were pregnancy or breastfeeding, any acute illness, history of congestive heart failure, atrial brillation (paroxysmal or persistent), myocardial infarction or stroke, neoplastic disease (past or current), thyroid disease, liver injury and chronic kidney disease (stage 3 or above) and smoking. Past smokers who stopped smoking at least 5 years before and accumulated less than 7 pack-years were permitted in the study. Asthma medications, except for biological treatment, were permitted, with oral corticosteroids at a daily dose ≤10 mg of prednisolone equivalent unless the doses were unchanged during the preceding 3 months.
Half of the subjects (N=20) were diagnosed with eosinophilic asthma de ned as having at least 1% of eosinophils in bronchoalveolar lavage uid (BALF), thus comprised the EA group. The remaining individuals formed the non-eosinophilic asthma (non-EA) group (N=20).

Spirometry studies
Spirometry with bronchial reversibility test using 400 mg of albuterol and body plethysmography after bronchodilator with assessment of the residual volume (RV) and total lung capacity (TLC) were carried out to assess differences in lung function between eosinophilic and non-EA patients using a Jaeger MasterLab spirometer (Jaeger-Toennies GmbH; Hochberg, Germany).
Microarrays were scanned with a InnoScan 900 Microarray Scanner and hybridization signals were detected using Mapix software (v.6.0.1; Innopsys, Carbonne, France). The software failed to perform microarray gridding correctly in some cases. Since gridding is crucial for image analysis as an error occurring during this process may be propagated to further stages 7 , a custom gridding algorithm to avoid grid displacement was developed. The coordinates of positioning markers to establish sub-grids were found using an MSER 8 -based method. Computed results were veri ed manually and corrected if needed, with overall gridding performance exceeding that of the commercial software, ultimately resulting in expression pro les of 33519 gene products. Microarray data preprocessing (background correction and normalization) was done by limma R software using the established methods 9,10 . Genomic data were analyzed with Bioconductor v.3.7. software of the R environment v.3.5.0.

Differential gene expression and co-expression identi cation and analysis.
A list of differentially expressed genes (DEGs) was established using limma package to assess for possible confounding with differential co-expression 11 ; false discovery rate (FDR) was assessed with Benjamini and Hochberg procedure. A list of differentially co-expressed genes (DCGs) was established using CoXpress 12 R plugin. In this method, groups of genes which expression levels are highly correlated in one set of experiments (i.e. EA), but not signi cantly correlated in a second set of experiments (i.e. non-EA) are found. A summary statistic of gene expression correlation matrix is made for each sample group and compared to distributions of n=1000 similar summary statistics generated by selecting groups of random genes of the same size from the same data set. The summary statistic used is the t-statistic taken from a test of the unique, pairwise, Pearson correlation coe cients against zero. The proportions of random statistics above the observed statistic are used to create a "probability of randomness" for each group in the two sample groups. A group of genes which distribution of pairwise correlations is found to be random in non-EA (p ≥0.05), and non-random in EA (p <0.05), is said to be differentially co-expressed between the two groups.
1.9 Gene set enrichment, creation of interaction networks and their topological robustness analysis.
Associations between DCGs, airway in ammation and remodeling were investigated using the following work ow: 1. Gene set enrichment analysis (GSEA) was performed to establish histone modi cations and cell types matching DCGs' expression pattern as well as associated biological processes, molecular functions and signaling pathways. The analysis was performed using Enrichr tool 13 using curated datasets (Tab. 1).
2. Based on curated resources ( 1), we established (1) lung-speci c binary protein-protein interactions of DCGs, (2) DCGs' upstream regulators and (3) DCG-microRNA-transcription factor interactions. Semiautomated analytics platforms Cytoscape 3.7.1 and NetworkAnalyst 3.0 were used to construct interaction networks as undirected graphs, where individual molecular entities formed the nodes and their interactions the edges.
3. Topological robustness analysis of the resulting graphs was performed by the previously described principle of topological attack on biological network 14 . Betweenness centrality was used to identify nodes with highest network-disruptive potential as possible targets for pharmacological interventions.
4. Betweenness centrality of the nodes was plotted against MEDLINE literature coverage using the keyword "asthma" to identify already established and candidate molecular entities of potential relevance. GenClip 3.0 sofware was used for automated literature mining.
5. Both individual DCGs and their rst-degree interactors were investigated in an attempt at mechanistic explanation of their role in the disease by literature mining using the keyword "asthma" and either of allergy and airway in ammation-related or airway remodeling-related keywords: (1) allergy and airway in ammation: Th1, Th2, Th17, dendritic cell, antigen-presenting, granulocyte, neutrophil, eosinophil, basophil, mastocyte, mast cell, monocyte, macrophage, lymphocyte, in ammation, thromboxane, prostaglandin, leukotriene, oxidative stress, allergy, allergic, atopy (1)  Demographic, clinical and laboratory characteristics of the study subjects are provided in Tables 2-4. EA and non-EA patients were comparable regarding most demographic and clinical features except for BMI, which was higher in the latter.
While the studied groups did not differ signi cantly in regard to lung function parameters, EA patients had more prominent absolute indices of air ow limitation (lower forced expiratory volume in 1 second; FEV 1 ) and more favorable obstruction reversibility (higher DFEV 1 ), while non-EA patients had lower total lung capacity (TLC), possibly due to overrepresentation of obese subjects in this group. HRCT indices of airway remodeling were signi cantly more prominent in the EA group, consistent with the differences in lung function, suggesting computed tomography as superior to spirometry in detecting advancement of airway remodeling.
2.2 Identi cation and study of differentially expressed genes.
After background correction, quantile normalization and ltration for coe cients of gene expression variation between 0.3 and 10, 14823 annotated gene products were included in subsequent analyses. Expression levels of 20 genes differed signi cantly between EA and non-EA by unadjusted p-value.
According to automated literature mining, at least one of them was described in the context of asthma alone in one paper, three in the context of in ammation-related keywords over 8 papers and 9 in the context of remodeling-related keywords over 80 papers (Tab. 5, Fig. 1).
After adjustment for false discovery rate, no gene's expression level was found to signi cantly differ between EA and non-EA, concluding DEG analysis and indicating unconfounded pure differential coexpression discovery 11 .
2.3 Identi cation and study of differentially co-expressed genes.
Hierarchical cluster analysis revealed a total of 23 groups of DCGs. A single group consisting of 32 genes with the highest mean, pairwise difference between the correlation matrices of gene expressions in the two asthma groups was chosen for subsequent analyses (Tab. 6).
Only 5 of the 32 genes co-occurred with asthma over 9 papers, 6 with allergy/in ammation-related keywords over 98 papers and 15 with remodeling-related keywords over 133 papers (Fig. 2). Abundance of signi cant correlations between DCGs' expression levels within EA group may indicate its pathogenetic consistency, contrary to less correlated DCGs in the non-EA group (Fig. 3). The resulting list of candidate genes as well as their interactors and regulatory networks ( Fig. 4-6) which may differentiate airway in ammation and remodeling mechanisms speci c to EA and non-EA are discussed.
None of the DCGs were signi cantly enriched for ontology terms at adjusted p-value (false discovery rate) < 0.05, expectedly, given the small gene set size. The results were ranked based on combined p-value and zscore.
2.4.1 Cell type-speci c histone modi cations related to expression of DCGs.
Fifteen of the 32 DCGs (all up-regulated in the EA) matched differential expression pattern related to histone modi cations H3K9ac and H3K27me3 of human lung broblast hg19, suggesting them as the likely source of the gene expression signal (Tab. 7). H3K9ac and H3K27me3 are the two histone modi cations of dominant and opposite role in regulation of neuroectodermal differentiation (H3K9ac) or pluripotency (H3K27me3) 15 . With H3K27me3 previously linked to epithelial-mesenchymal transition (EMT) and extracellular matrix degradation upon treatment of cell culture with TNF-a and TGF-b 16 , up-regulation of the genes related to both histone modi cations in our study could indicate their role in ongoing EMT in EA, consistent with the abundance of eosinophil-derived TGF-b1 in this disease phenotype 17 .

Gene ontology associations.
Molecular functions and biological processes associated with DCGs by gene ontology are listed in Tab. 8-9. They include phenomena signi cant in the pathogenesis of asthma and airway remodeling, such as Tcell receptor binding, MHC II protein complex binding, immunoglobulin secretion, myeloid cell apoptosis, regulation of cellular response to growth factors (e.g. TGF-b) and metalloendopeptidase inhibitory activity, as further discussed in section 3.5.
These kinase-dependent DCGs might indicated differential activity of the related pathways and have pathogenetic role in EA. The ligand of TGFBR2, TGF-b, has well-established role in asthma and airway remodeling, as further discussed in section 3.5, while the receptor itself was found to be involved in T-cell differentiation. The role of HIPK1 was to date only described in fetal angiogenesis activated by TGF-β-TAK1 pathway 18 , apoptosis induced by TNF-a and prevention of MAP3K5-JNK activation. Our study links HIPK1 with eosinophilic asthma through the pattern of DCGs, a suggestion reinforced by association between MAP3K5 and atopy 19 . IKBKE is a known target of the nuclear factor κB (NFκB), a transcription factor involved in chronic in ammation of established role in asthma pathogenesis, including remodeling.
2.5 Literature-based assessment of individual DCGs for mechanistic role in EA and bronchial remodeling.
Individual DCGs are discussed below in the following order: gene's function, possible involvement in the disease and interaction with other DCGs from our study (shown in bold).
2.5.1 Sodium/potassium-transporting ATPase subunit beta-1 (ATP1B1) The sodium-potassium pump is a heterodimer composed of the catalytic subunit alpha and non-catalytic subunit beta. The latter may be expressed in several cell types and conditions involved in asthma pathogenesis: (1) In the epithelia, the beta-1 subunit may be related to epithelial sheathing as it contributes to formation and stabilization of intercellular junctions. It regulates the number of pumps transported to plasma membrane through assembly of alpha/beta heterodimers, which can act as positive or negative regulator of intracellular adhesion, depending on its proportion to FXYD5, a membrane protein involved in chemokine up-regulation and loosening of cell adhesion through down-regulation of E-cadherin 20 .
(2) In cytomegalovirus (CMV) infection, previously linked to asthma pathogenesis via promotion of Th2 response 2122 , subunit beta 1 was found to co-localize with viral UL136 protein suggesting involvement in cell-to-cell spread of the infection 23 . It is also upregulated in human neutrophils during respiratory syncytial virus (RSV) infection, a trigger of asthma exacerbations.
(3) In lymphocytes, subunit beta 1 is known to be up-regulated in response to phytohemagglutinin, suggesting possible role as marker of their activation status 24 and AT1B1 was among four genes upregulated in human lymphocytes under Th17-polarizing conditions.
(4) Decreased activity the pump elevates intracellular free calcium concentration, a phenomenon linked to the activation of in ammatory cells and airway smooth muscles (ASM) resulting in increased proliferation, spreading, and eosinophil-attracting eotaxin-1 release in experimental conditions 25262728 .
(5) In the goblet cells, IL-4 produced by Th2-lymphocytes causes up-regulation and translocation of ATP1B1 from basal to apical aspect of the cell, where it colocalizes with H+/K+-ATPase, ATP12A, being required for its function implied in mucous secretion in asthma and cystic brosis 29 .
(6) ATP1B1 is modulated by PXK, which is thought to participate in regulation of electrical excitability and synaptic transmission by, in part, epidermal growth factor (EGF) receptor ligand-induced internalization 30 , both processes being postulated to play a role in airway hyperreactivity in asthma. ATP1B1, CLC, FADS6, FBN3, RECK and SLC19A1 share a common transcription factor Sp1 313233 involved in multiple cellular processes ranging from cell proliferation and differentiation to immune response, their speci c function being highly dependent on post-translational modi cation. Sp1 expression and DNA binding activity are increased in CMV infection 34 , a positive-feedback loop given ATP1B1's role in cell-tocell virus transmission. Furthermore, Sp1 has an established role in asthma pathogenesis as a transcription factor for (1) lipooxygenase ALOX5 35 , which promoter region's polymorphism affects asthma control, and (2) vascular endothelial growth factor (VEGF), hypersecreted from ASM in asthma and contributing to remodeling. Of note, both ATP1B1 and RECK are implied to take part in airway remodeling in mustard lung 36 , adding to evidence of Sp1-regulated genes in airway pathology. ATP1B1 was one of the genes included in transcriptomic cluster TAC3 of patients with moderate to high sputum eosinophilia in the U-BIOPRED cohort of Th2-dependent asthma 37 , with our study being the second to indicate the link.
Of note, as ATP1B1 is involved in kidney proximal tubule bicarbonate reclamation, aldosterone-regulated sodium reabsorption and kidney stone formation, common regulation of ATP1B1 in both bronchi and renal tubule may underlie apparent higher incidence of nephrolithiasis and chronic kidney disease in asthma patients 3839 .

Charcot-Leyden crystal protein (CLC)
Charcot-Leyden crystal protein or galectin-10 (Gal-10) is an atypical galectin preferentially binding bmannosides. It takes part in sequestration and vesicular transport of eosinophil granule cationic ribonucleases: eosinophil-derived neurotoxin and eosinophil cationic protein 40 and was found to poses lysolecithin acylhydrolase (lysophospholipase; phospholipase B) activity. Released upon degranulation, it is ubiquitous in sputum of EA patients and considered a marker of eosinophilic in ammation. It is thought to regulate immune responses through the recognition of cell-surface glycans and may possess IgEbinding capability. Furthermore, through lysophospholipase activity, it may hydrolyze 1) phospholipids, producing arachidonic acid for synthesis of eicosanoids of important role in asthma pathogenesis and 2) surfactant phosphatidylcholine, leading to surfactant dysfunction and small airway closure 41 .
Recently, the existence of regulatory CD16-high eosinophils with distinct suppressor function on T-cell proliferation has been suggested, with their function likely mediated by Gal-10 in immune synapses between eosinophils and lymphocytes 42 . Furthermore, Gal-10 was recently found to be expressed in human CD4(+)CD25(+)Foxp3(+) T regulatory cells (Treg), while nearly absent in resting and other activated CD4(+) T cells. In Treg cells this lectin is essential to limit proliferation and suppressive function, similarly to the above CD16-high eosinophils 43 .
CLC shares transcription factor Sp1 with ATP1B1, RECK, FBN3, FADS6 and SLC19A1, which points to a possible involvement in a common biological process.
Consequently, CLC was included in transcriptomic cluster TAC1 of patients with high sputum eosinophilia in the abovementioned U-BIOPRED asthma study 37 and associated with asthma in an epigenome-wide association study (EWAS) 44 .
Since galectins can frequently act as both positive and negative modulators of the same processes, possibly depending on isoform plasticity, posttranslational modi cations, localization or cofactors 45 , further research is needed to assess if galectin-10 can be in uenced to act as an anti-in ammatory factor, similar to galactin-9 (Gal-9). This galectin, released by various cell types plays a substantial role in control of effector cells and Th1/Th2 balance. It is known to be upregulated upon CMV infection through induction of IFN-b 46 and released upon calcium-mediated exocytosis 47 . It is able to reduce Th2-associated airway in ammation and airway hyperresponsiveness through binding with T cell immunoglobulin mucin domain 3 (Tim-3), persistently expressed on functional T-cells during chronic viral infections. By inducing maturation of monocyte-derived dendritic cells it promotes Th1-associated immune response. It also triggers apoptosis of mature Th1 cells by aforementioned Tim-3 receptor, preventing degranulation of mast cells by forming complex with IgE 48 and inhibits interaction-dependent migration of in ammatory cells via CD44 pathway [49][50][51][52] . These properties elicit galectin-9 and possibly galectin-10 as potential candidates for development of novel protein drugs for the treatment of asthma. Co-expression of the two genes involved in limitation of Th2-and induction of Th17-differentiation, CLC and PSG2, in EA patients suggests both a counterweight mechanism to limit Th2-response and a possible role behind heterogeneity of Th2-high asthma, recently proposed to be divided into IL-5-high/IL-17F-high asthma (with mixed granulocytic in ltration) and IL-4/IL-13-high asthma (with eosinophilic in ltration alone) 56 .

EPS8 like 1 (EPS8L1)
EPS8L1 encodes a protein related to the epidermal growth factor receptor kinase substrate 8 or Eps8 involved in actin remodeling. EPS8L1 itself is known to play a role in T-cell receptor binding, membrane ru ing and remodeling of actin cytoskeleton through F-actin organization by stimulating guanine exchange of SOS1, thus taking part in regulation of cell locomotion. In a murine knockout models, EPS8L1 was required for EGF-dependent membrane ru ing 57 and Eps8 was required for maintaining front-to-back polarity of dendritic cell, both features required for cell migration 58 . Furthermore, Eps8 takes part in Cdc42-IRSp53-Eps8 pathway involved in formation of lopodia and cell migration, with EPS8L1 being the only Eps8-related protein associating with IRSp53 and thus possibly involved in the process 59 .
No direct link between EPS8L1 an asthma was made to date.
Cell migration is a regulated by Rho GTPases, of which RhoBTB2 is an example, although no associations between Eps8 or EPS8L1 exist to date.

Rho related BTB domain containing 2 (RHOBTB2)
RHOBTB2 is an atypical member of the Rho family of small GTPases which control cell migration, invasion and cycle and is required for expression of CXC motif ligand 14 (CXCL14). CXCL14 is a chemoattractant that controls dendritic cell activation, leukocyte migration and angiogenesis and an autocrine growth factor for broblasts facilitating their migration 60 and which expression is lost through unknown mechanisms in a wide range of epithelial cancers 61 ; e.g., loss of RhoBTB2 expression correlates with downregulation and loss of CXCL14 secretion by head and neck squamous cell carcinoma cell lines, whereas reintroduction of RhoBTB2 restores CXCL14 secretion 62 .
In a genome-wide association study conducted in Australia, CXCL14 was one of the six most-associated loci with the asthma susceptibility 63 . Moreover, effects of exposure to tobacco smoke pollution are thought to be mediated through the CXCL14. The experimental data suggests differential effects of occasional 64 and chronic 65 tobacco exposure on CXCL14 expression, as well as differential effects of tobacco smoking on airway eosinophilia (elevated in chronic smokers, but attenuated by occasional exposure in non-smokers).
While CXCL14 receptor is not yet identi ed, migration of antigen-presenting cells (APCs) towards CXCL13, the ligand of CXCR5, is signi cantly potentiated in the presence of CXCL14, possibly contributing to mucosal recruitment of in ammatory cells 66 , including APCs and Th17 lymphocytes expressing CXCR5 as well. Regulation by RhoBTB2 suggests its role in disease mediated through CXCL14 and resulting migration of leukocytes, dendritic cells and possibly epithelial-mesenchymal transition of broblasts.
Furthermore, RhoBTB2 is a substrate adaptor for Cul3-based ubiquitin ligase complex known to play a role in a Keap1-guided degradation of antioxidative Nrf transcription factor implicated in asthma pathogenesis 67 . As Keap1 is believed to compete with other BTB proteins for culling-3 binding, downregulation of RhoBTB2 may leave E3 ubiquitin ligase open for adaptor Keap1-BTB-mediated ubiquitination and degradation of antioxidative Nrf2, thus suppressing its transcriptional activity and promoting oxidative stress. Recent studies indicate that RhoBTB2 is involved in the Hippo signaling pathway through LKB1 regulation, with loss of RhoBTB2 leading to ubiquitination and loss of LKB1 as well as increased YAP activity of role in asthma pathogenesis. Available GEO dataset indicates correlation between RhoBTB2 expression and asthma exacerbation, with signi cant difference in expression between pools of peripheral blood mononuclear cells from exacerbated and convalescent asthma patients (GDS3615 68 ). Furthermore, RhoBTB2 was one of the genes which expression changed signi cantly after bronchial thermoplasty 69 .
RhoBTB2's GTPase activity as well as interaction with Cul3 requires Hsp90, a chaperone regulated by CCT7 70 . Its role in membrane ru ing makes a possible link with aforementioned EPS8L1.
2.5.6 Ras association (RalGDS/AF-6) and pleckstrin homology domains (RAPH1) RAPH1 is a member of Mig10/Rap1-interacting family of adaptor proteins regulating actin dynamics. As a mediator of localized membrane signaling implicated in the regulation of lamellipodial dynamics it regulates cell migration, and its internalization negatively regulates cell adhesion, possibly contributing to epithelial sheathing. Phosphorylation of RAPH1 by C-abl oncogene takes part in its coordination of actin remodeling in response to external stimuli 61 . As a target gene of regulatory miR-203, it facilitates keratinocyte migration and wound healing 71 . While no direct role of RAPH1 in airways was reported, miR-203 is known to inhibit proliferation in ASM cells through attenuation of c-Abl expression and phosphorylation by ERK1/2, a pathway involved in epithelial-mesenchymal transition in airways, of which RAPH1 is a phosphorylation target. Furthermore, both ERK1/2 and RAPH1 are affected by histone demethylase KDM2A, implicated in lung carcinogenesis 72 . Because of its place at the intersection of KDM2A, miR-203, c-Abl and ERK1/2, RAPH1 upregulation in EA patients suggests its role in ASM proliferation, as well as epithelial sheathing and epithelial-mesenchymal transition.
2.5.7 Signal recognition particle receptor beta subunit (SRPRB) SRPRB is one of two subunits of the signal recognition particle receptor (SRP) required for co-translational targeting of secretory and membrane proteins. With the alpha subunit targeting SRP-ribosome-nascent polypeptide complexes to translocon, the beta subunit is a transmembrane GTPase which anchors the alpha subunit to the endoplasmic reticulum (ER).
ILK1, an integrin-linked kinase responding to signals from the extracellular matrix, including injury, and regulating basal stem cell activity and ASM differentiation 74 links SRPRB to CCT7 through a nity capture protein-protein interaction of unknown biological signi cance.

Ankyrin repeat domain 26 pseudogene 1 (ANKRD26P1)
ANKRD26P1 is a pseudogene of little known function. It interacts with MAGEA-6 75 , a protein encoding epitopes presented with HLA-DRB1*0401. Presentation of these epitopes recognized by CD4(+) T-cells from patients with melanoma or renal cell carcinoma enhances T-cell response against cells expressing both proteins in HLA-DRB1*0401 patients 76 . The same haplotype was previously associated with severity of asthma, incidence of rheumatoid arthritis and eosinophilic granulomatosis with polyangiitis, an autoimmune disease resulting in eosinophilic tissue in ltrates and EA phenotype in its clinical presentation 77,78 .
MAGEA-6 is up-regulated upon exposure of cell culture to macrophage migration inhibitory factor (MIF), a pro-in ammatory factor known to play a role in asthma and chronic obstructive pulmonary disease 79 .
These ndings suggest that genetic factors such as human leukocyte antigens may be associated with susceptibility to EA and that it may be mediated through T-cell response against cells expressing ANKRD26P1 and MAGEA-6, possibly in response to MIF and resulting in T-cell production of IL-4 and IL-13 that drive eosinophilic in ammation.

Chaperonin containing TCP1 subunit 7 (CCT7)
CCT7 is a molecular chaperon, member of chaperonin containing TCP1 complex (CCT; TRiC) assisting folding of proteins like actin and tubulin and regulating Hsp90 chaperone.
Its expression increases in brotic wound healing and is essential for the accumulation of α-smooth muscle actin (α-SMA) in broblasts and cell differentiation to myo broblasts 80 . CCT7 was recently found to coimmunoprecipitate with thromboxane A2 receptor as well as β2-adrenergic receptor from human HEK 293 cells and its depletion resulted in reduced cell surface expression of both receptors 81 ; a fact signi cant given the momentous role of both G-protein coupled receptors in asthma pathogenesis. Two-hybrid screening revealed interaction between CCT7 and ANXA1, which de ciency causes spontaneous airway in ammation and hyperresponsiveness in murine studies 82 .
Notable is the CCT7's central position in the PPI network of multiple DCGs and their interactors of high asthma literature coverage (Fig. 4) and involvement in epithelial-mesenchymal transition, implying regulatory function of the CCT/TRiC complex in asthma pathogenesis, involving e.g. Hsp90-dependent GTPase activity of RhoBTB2.
TGF-b is known to induce DIO3 expression in various cell types, including lung broblasts. Of note, deiodinase type III is a selenoprotein while selenium de ciency has been suggested to take part in in ammatory diseases including asthma 83 through contribution to oxidative stress, a mechanism thought to underlie asthma exacerbations in hyperthyroidism. As oxidative stress is implied to induce DIO3 expression, asthma exacerbations in selenium de ciency are aggravated through impaired DIO3 function and increased tissue exposure to active thyroid hormones. No direct link between DIO3 and asthma has been made to date, however, its upregulation might be a response to oxidative stress, activation of fetal genes due to tissue repair (along with FBN3 and GPI, below) or adult overexpression of genes hypermethylated prenatally, contributing to the hypothesis of prenatal epigenetic changes as predisposing to asthma development 52 . Promotor of the above-mentioned ATP1B1 has sequences binding thyroid hormones (TRE), forming a link to upregulation of both ATP1B1 and DIO3 in our study.
2.5.11 Epidermal growth factor receptor pathway substrate 15 (EPS15) EPS15 is a protein involved in the clathrin-dependent internalization of ligand-inducible tyrosine kinase receptor (RTK) types (including EGFR, TGF-b receptors and integrin b-1) as well as cell adhesion molecules (integrins, E-cadherin), receptors relevant for bronchoconstriction (β2-adrenergic and M3-muscarinic receptors) and antigen-presenting major histocompatibility class II proteins.
EPS15 is also necessary for non-clathrin-dependent endocytosis of EGFR and TGF-b, with knockout models resulting in attenuated TGF-β-induced Smad2 phosphorylation 84 of role in broblast to myo broblast transition 85 .
EPS15 forms PARK2-EPS15-EGFR complex with parkin, an E3 ubiquitin ligase in uenced by IL-13 and known to enhance airway mitochondrial DNA release contributing to in ammation 86 . Since IL13-stimulated parkin aids mono-ubiquitination of EPS15 thus preventing RTKs internalization and potentializing the effects of TGF-b, it forms a link between the two pro-brotic cytokines as well as mtDNA-mediated in ammation.
EPS15L1, a protein of suggested role in the process of receptor endocytosis, coimmunoprecipitates with EPS15 and is one of the genes associated with asthma in an EWAS study 44  FADS6 is regulated by hsa-miR-331-3p, a micro-RNA post-transcriptional regulator associated with lung function in asthma, with our study supplying evidence for the mechanism of such association being mediated through FADS6 94 . FADS6 promoter binds factors reported to in uence asthma risk: EZH2, as well as abovementioned Sp1, a transcription factor for ATP1B1, FBN3, RECK and SCL19A1.

Fibrillin 3 (FBN3)
Fibrillins are components of extracellular calcium-binding micro brils, which occur either in association with elastin or in elastin-free bundles and act as structural support important for extracellular matrix integrity. Like PSG2, brillins may be involved in regulation of TGF-b activity through association with latent TGF-b binding proteins. Physiological FBN3 expression was described only in fetal developing bronchi 9596 , but its mutations have been found in lung cancer 97 . It shares a common transcription factor Sp1 with ATP1B1, FADS6, RECK, SLC19A1 33 . Up-regulation of the FBN3 gene in the current gene set might both reinforce the hypothesis of fetal program execution during bronchial remodeling (together with GPI and DIO3), regulation of signaling through latent TGF-b and add to a weak association between lung cancer and asthma 98 .
GPI has a proven pathogenetic role in rheumatoid arthritis, acting as an autocrine synovial broblast proin ammatory cytokine resulting in proliferation, apoptosis inhibition and secretion of tumor-necrosis factor (TNF)-a and IL-1β 99 . Immunization of murine CD4(+) T-cells against GPI causes experimental model polyarthritis with production of proin ammatory TNF-a, IL-17 and IL-6 and anti-GPI autoantibodies are detectable in some human rheumatoid arthritis patients. GPI was found to be a direct regulator VEGFmediated angiogenesis in rheumatoid arthritis, with both proteins being up-regulated by HIF-1a, a transcription factor involved in airway remodeling. Though rheumatoid arthritis is classically considered as Th1-dependent disease, it seems to share genetic risk factors and elements of pathogenesis with asthma 100 . Further studies are needed to assess if anti-GPI is present in sputum of asthma patients similarly to anti-eosinophil peroxidase and anti-IgE autoantibodies and if VEGF-mediated airway remodeling is regulated by GPI in similar fashion to that seen in rheumatoid arthritis.
As an enzyme taking part in glycolysis, it may affect glucose metabolism in the airway epithelia, where hyperglycemia and impaired glucose transport induce cell proliferation 101 . Association with asthma in an EWAS 44 adds evidence of its pathogenetic role. Its role in fetal development along with DIO3 and FBN3 may be indicative of fetal genetic program execution in asthmatic airways.

Mitochondrial ribosomal protein L14 (MRPL14)
MRPL14 is a nuclear gene that encodes part of two intersubunit bridges in the assembled mitochondrial ribosome required for mitochondrial translation. Related to asthma in a genome-wide association study 102 and controlled by MYC transcription factor similarly to DGLUCY, it may be involved in asthma-speci c mitochondrial biogenesis and accumulation in ASM.

Olfactory receptor family 52 subfamily I member 1 (OR52I1)
OR52I1 is an olfactory receptor functioning by CaBP/Cam-dependent inhibition of calcium channel, a function which may be in uenced by CABP5 due to functional replaceability with calmodulin. While no link between OR52I1 and asthma was found to date, other olfactory receptors like OR1D2 and OR2AG1 are known to be expressed in ASM cells and their activation by amyl butyrate was shown to inhibit histamine-induced ASM contraction. OR52I1 and CABP5 may thus be involved in pathology of airway hyperreactivity in asthma of undetermined signi cance.

Protein phosphatase 2 regulatory subunit B beta(PPP2R3B)
PPP2R3B might modulate its substrate selectivity and catalytic activity and direct the localization of the catalytic enzyme to a particular subcellular compartment. Also known as PR48, it regulates cell cycle progression likely by controlling initiation of DNA replication 103 . It's downregulation in EA may indicate increased airway or in ammatory cell proliferation.

Reversion inducing cysteine rich protein with kazal motifs (RECK)
RECK is a membrane-anchored negative regulator of matrix metalloproteinase-9 (MMP)-9 able to suppress its secretion and enzymatic activity. MMP-9, a gelatinase, is known to be secreted by epithelial cells, neutrophils and eosinophils in response to allergen challenge or TNF-a and was previously linked to both eosinophilic in ltration into the bronchial wall and remodeling as a major MMP involved in asthma 104 . Its activity is increased in sputum of asthmatic patients and correlates with loss of FEV1 in response to allergens 105 . Furthermore, RECK down-regulation by oncogenic signals may facilitate tumor invasion and metastasis through regulation of MMP-2 and MT1-MMP, further linking asthma and carcinogenesis. It shares a common transcription factor Sp1 with ATP1B1, FADS6, FBN3 and SLC19A1 32 .

Scavenger receptor cysteine rich family member with 4 domains (SSCRB4D)
SSCRB4D belongs to the scavenger receptor cysteine-rich (SRCR) superfamily of highly conserved proteins involved in the development of the immune system and regulation of immune responses. SSCRB4D may be associated with asthma by being one of the target genes of BACH1, involved in the response to oxidative stress and previously associated with asthma in an EWAS study 44 .

Stomatin like 3 (STOML3)
STOML3 modulates mechanotransduction channels and acid-sensing ion channels (ASIC) proteins and potentiates PIEZO1 and PIEZO2 function by increasing their sensitivity to mechanical stimulations. Studies in animal models suggest that mechanotransduction may play a role in hyperresponsiveness of the airways in asthma, while PIEZO1 de ciency decreases cell adhesion and increases migration of lung epithelial cells 106 , a feature of asthmatic bronchial remodeling 107 .

Striatin-4 (STRN4)
STRN4 belongs to a family of proteins binding calmodulin in a calcium-dependent matter, which may function as scaffolding or signaling protein. It is known to co-localize with protein phosphatase 2A (PP2A) acting as a regulatory subunit of the STRIPAK kinase-phosphatase complex involved in multiple cellular processes. Animal studies suggest that STRIPAK might be involved in cell cycle control, cell adhesion, migration, epithelial integrity and epithelial-to-mesenchymal transition. STRN4, by inhibition of MAP4K4 kinase, is considered the key regulator inhibiting the Hippo pathway involved in asthma pathogenesis 108 . In T-cells, de ciency of MAP4K4 results in Th17 differentiation with IL-6 and IL-17 expression 109 . In a study of asthmatic patients, miRNAs expressed in the airways frequently targeted MAP kinases including MAP4K4 110 , suggesting their involvement in the disease.
STRN4 itself was also proposed as a part of mechanism to control dendritic spine morphology, adding to possible neuronal involvement in airway disease. STRN4 is one of the genes associated with asthma in an EWAS 44 .

Tetratricopeptide repeat domain 3 pseudogene 1 (TTC3P1)
TTC3P1 is a pseudogene resembling TTC3, a E3 ubiquitin-protein ligase that mediates ubiquitination of phosphorylated Akt (AKT1, AKT2 and AKT3) as part of negative feedback required to control Akt levels after pathway activation. TTC3 contributes to TGF-b1-induced epithelial-mesenchymal transition any myo broblast differentiation, forming a direct link with bronchial remodeling 112 . TTC3 may also play a role in neuronal differentiation inhibition via its interaction with CIT, associated with asthma in GWASdb database 102 , a fact of potential signi cance for the neuronal involvement in the disease 113 . As an apparently transcribed pseudogene, TTC3P1 may regulate TTC3 function. Membership in a E3 ligase complex associated with CUL3 links TTC3 to RhoBTB2 while involvement in Akt pathway makes an association with ASB3.

Calcium binding protein 5 (CABP5)
CABP5 product inhibits calcium-dependent inactivation of L-type calcium channel, thus increasing cell excitability. While primarily involved in the transmission of light sensation in the retina, it was reported in Tcells including Th1, Th17, mucosa-associated invariant T (MAIT) cells and B-cells, being differentially expressed in allergic asthma patients 114 ). Together with the involvement in stimulating neurite outgrowth and vesicle exo-as well as endocytosis in PC12 cells, these facts suggest that CABP5 plays a role in immune reactions, neuronal growth (a feature of airway remodeling) and vesicle handling. Able to functionally replace calmodulin, it may potentially interact with striatin in clathrin-dependent membrane receptor endocytosis, thus relating CABP5 to STRN4 and EPS15 115 .

Ribosomal protein S13 (RPS13)
RPS13 is a component of the 40S ribosomal subunit essential for eukaryotic protein synthesis and regulated through redundant mechanisms, hence described as a reference housekeeping gene. However, it is known to take part in G1 to S cell phase transition and was featured in cluster 1 of asthma-linked modules obtained from the combination of differential gene expression and GWAS study of the U-BIOPRED project 116 . Furthermore, it forms TNF-a/NF-Kappa B signaling complex of established role in asthma 67 .
These data provide arguments for viewing RPS13 as an active player in asthma pathogenesis and reevaluate it as a reference housekeeping gene of seemingly stable expression across biological conditions. A nity capture studies indicate interactions of RPS13, CCT7 and EPS15 with VCAM1, ITGA4, FN1 proteins involved in cell-cell and cell-extracellular matrix adhesion, with VCAM and FN1 both binding to ITGA4, as subunit of integrin a4b1. While the nature of the interaction between proteins involved in extracellular matrix and endocytosis, regulation of chaperones and ribosome structure is unknown.

Ankyrin Repeat And SOCS Box Containing 3 (ASB3)
ASB3 is part of ASB gene family involved in Erk1/2 and PI3K/Akt signal transduction pathways by regulation of MAP kinase and Akt phosphorylation, both implicated in smooth muscle proliferation in asthma 117 .
It may regulate lymphocyte differentiation through TNF-a receptor (TNFR) 2 ubiquitination by recruitment of E3 ubiquitin ligase adaptors: elongins-B/C 118 . While TNFR2 initiates immune modulation and tissue regeneration, signaling through TNFR1 triggers pro-in ammatory pathways 119 . Loss of TNFR2 signaling can impair expansion and stability of Tregs and decrease their sensitivity to IL-2. Moreover, through reciprocal PI3K/Akt pathway activation and phosphorylation of STAT5, ASB3 impairs Th17 differentiation and may impair IL-12 signaling through JAK2/STAT4 to favor Th2 rather than Th1 differentiation 120 , resulting in predominantly eosinophilic, rather than neutrophilic in ammation.
Recent GWAS study has associated three SNPs in a region of chromosome 2 near ASB3 and SOCS with degree of bronchodilation following inhalations of albuterol in asthma patients 121 , with our study indicating ASB3 as an agent in the disease pathogenesis.

Solute carrier family 19 member 1 (SLC19A1)
SLC19A1, also referred to as reduced folate carrier protein, is a bidirectional membrane anion transporter enabling both folate and anti-folate medications, e.g. methotrexate (MTX) in ux into the cell. Folate is required for synthesis of purines and thymine, homocysteine-methionine metabolism, both homocysteinedependent and independent production of nitric oxide and reactive oxygen species (ROS) as well as DNA methylation, exerting both anti-in ammatory and pro-in ammatory properties, possibly depending on chronicity of an in ammatory disease 122 .
Folate de ciency was previously linked to asthma severity 123 as well as risk of developing the disease in children of folate-de cient mothers 124 . It affects broblast expression of genes related to cytoskeleton remodeling, extracellular matrix and signaling through Wnt pathway (DKK1, WISP1 and WNT5A) 125 and may affect collagen metabolism in concomitant ascorbic acid de ciency. Folate de ciency increases the CD4/CD8 ratio and murine models indicate its role in maintaining CD4(+)Foxp3(+) regulatory T-cells, decreased in human asthma patients 126 .
Furthermore, as a recently discovered major transporter of immunoreactive cyclic dinucleotides (CDN) of autologous or microbial origin 127 , SLC19A1 may take part in dysregulation of the cGAS-STING pathway which overactivation can drive in ammatory lung diseases, including asthma 128 .
A genetic risk factor of atopy, an SNP rs12483377-A is mapped to SLC19A1 and COL18A1, both functionally implied in the process 19 .
Up-regulation of SLC19A1 may be driven by transcription factors linked to asthma (CEBPB, USF1, SP1), Th2-differentiation (CEBPB) 122 and asthma-speci c mitochondrial biogenesis (NRF-1). Furthermore, it is up-regulated by vitamin D 129 , which de ciency was previously linked to asthma control. Differential coexpression in EA may indicate both activity of the above transcription factors, cell proliferation status or impaired reaction to folate de ciency, a condition under which SLC19A1 is considered cytotoxic, since it may aggravate folate deprivation by expulsion of folate monoglutamates 130 . These facts indicate that SLC19A1 may have pathogenetic role in asthma, being affected by both transcription factors involved in asthma and nutritional status linked to asthma risk and severity.
We hypothesize that ill-regulated overexpression of SLC19A1 by activity of transcription factors to asthma pathogenesis in a bystander effect may sensitize cells to pro-in ammatory effects of self-and microbial DNA or aggravate cell folate deprivation, activating cGAS-STING pathway or creating a vicious cycle of abnormal DNA methylation and exhibition of folate's pro-in ammatory properties.
Co-expression of SLC19A1 in EA patients may have therapeutic implications, since MTX is consideration as steroid-sparing medication for severe asthma patients and its transport may affect treatment outcomes. MAEA is an ubiquitin ligase and part of CTLH complex of signi cance in TNF-a-mediated apoptosis, as well as erythropoiesis and macrophage maturation. An EWAS found association between trans-CpG site in MAEA and IL1RL1a, an isoform of IL1RL1 associated with multiple immune reactions in the lung, including Th2-response 131 . MAEA methylation positively correlated to birth weight in an observational study 132 , suggesting possible mechanistic association between obesity, MAEA and IL1RL1a activity in asthma. As an apoptosis-related gene it is up-regulated upon exposure to cigarette smoke extract, possibly due to circulating TNF-a, placental growth factor or loss of VEGF signaling, as well as RSV exposure, along with MMP-9 133 .
DGLUCY can be linked to asthma through (1) impact on E-cadherin expression (involved in to epithelial integrity in asthma), possibly through inhibition of ERK and P90RSK phosphorylation 134 , (2) its suppression by miRNA-199b (overexpressed during bacterial but not viral infections) 135 , (3) possibly related differential expression in neutrophilia 136 , (4) strong differential expression implied in seasonal remodeling of the immune system (reinforcing its role in immune system function and possibly regulation) 137 .
Expression of DGLUCY is in uenced by signaling through oestrogen receptor alpha 137

Asthma ignorome:
We couldn't nd relevant associations with asthma, in ammation, airway remodeling or other DCGs for GOLGA2P3Y, RP11-321E2.2 and RP3-473B4.3. They will become part of asthma ignorome 141 of currently unknown role until proven otherwise by future studies.
2.6 Regulatory networks of differentially co-expressed genes and inferred pathways.
Regulatory networks of DCGs are presented in Figures 5-6. Associated transcription factors were plotted against betweenness centrality of combined regulatory network to identify highly connected nodes with low literature representation and thus presenting research opportunities (Fig. 7). An enrichment analysis of the transcription factors revealed multiple signaling pathways listed in Tab. 11. Exceptional enrichment of SMAD2/3 nuclear signaling regulation corresponds to pro brotic effects of TGF-β/Smad2/3 pathway involved in broblast to myo broblast transition 85 . Activator protein-1 pathway is necessary for transcription of Th2-pro le cytokine IL-4 and indicated as therapeutic target in asthma 142 c. The E2F/Rb and Wnt/b-catenin pathways are both implied in airway smooth muscle hypertrophy 143,144 . This concordance with current literature validates our results and suggests other less studied pathways as important for future studies.

Summary.
Systems biology approach allowed us to place 29 of the 32 differentially co-expressed genes in a pathogenetic network with molecular entities previously described as involved in asthma pathogenesis, thus indicating novel candidates for asthma research and highly-connected network disruptors of potential therapeutic applications.
Genes differentially co-expressed in our study bear signi cance in airway viral infection (