Reticular Basement Membrane Thickness Is Associated with Growth- and Fibrosis-Promoting Airway Transcriptome Profile-Study in Asthma Patients

Airway remodeling in asthma is characterized by reticular basement membrane (RBM) thickening, likely related to epithelial structural and functional changes. Gene expression profiling of the airway epithelium might identify genes involved in bronchial structural alterations. We analyzed bronchial wall geometry (computed tomography (CT)), RBM thickness (histology), and the bronchial epithelium transcriptome profile (gene expression array) in moderate to severe persistent (n = 21) vs. no persistent (n = 19) airflow limitation asthmatics. RBM thickness was similar in the two studied subgroups. Among the genes associated with increased RBM thickness, the most essential were those engaged in cell activation, proliferation, and growth (e.g., CDK20, TACC2, ORC5, and NEK5) and inhibiting apoptosis (e.g., higher mRNA expression of RFN34, BIRC3, NAA16, and lower of RNF13, MRPL37, CACNA1G). Additionally, RBM thickness correlated with the expression of genes encoding extracellular matrix (ECM) components (LAMA3, USH2A), involved in ECM remodeling (LTBP1), neovascularization (FGD5, HPRT1), nerve functioning (TPH1, PCDHGC4), oxidative stress adaptation (RIT1, HSP90AB1), epigenetic modifications (OLMALINC, DNMT3A), and the innate immune response (STAP1, OAS2). Cluster analysis revealed that genes linked with RBM thickness were also related to thicker bronchial walls in CT. Our study suggests that the pro-fibrotic profile in the airway epithelial cell transcriptome is associated with a thicker RBM, and thus, may contribute to asthma airway remodeling.


Introduction
Asthma is an inflammatory disease of the airways that may result from exposure to inhaled allergens or other environmental agents, rendering genetically susceptible individuals to aberrant activation and prolonged immune responses [1]. Airway remodeling is a pathological feature of persistent asthma and contributes to the clinical manifestations of the disease. It refers to the bronchial wall's structural changes caused by chronic inflammation, repeated cycles of injury, and repair [2]. In asthma, airway remodeling is characterized

Patient Description and Basic Laboratory Tests
We analyzed 40 non-smoking moderate to severe asthma patients aged 30-65 years, with at least a 3-year history of confirmed asthma and no exacerbation within at least six months before the enrollment. About half of the patients (n = 21, 52.5%) had persistent airflow limitation, while the remaining (n = 19, 47.5%) had normal spirometry before or after bronchodilator. Persistent airflow limitation was defined as a forced expiratory volume in 1 s (FEV 1 )/vital capacity (VC) index below 0.7 or FEV 1 lower than 0.8 of predicted value after bronchodilator.
Both subgroups did not differ in the demographic parameters or asthma severity and duration. Patients with persistent airflow limitation were characterized by higher blood and bronchoalveolar lavage (BAL) eosinophilia (Table 2). Intriguingly, they also had lower interleukin (IL)-6 concentrations in BAL.    Table 3 presents cross-sectional geometry parameters of the right upper lobe apical segmental bronchus (RB1) and the right lower lobe basal posterior bronchus (RB10).   Table 3 presents cross-sectional geometry parameters of the right upper lobe apical segmental bronchus (RB1) and the right lower lobe basal posterior bronchus (RB10). As expected, corresponding airway cross-sectional geometry parameters measured in RB1 and RB10 correlated well with each other; for wall area ratio (WAR) and wall thickness ratio (WTR), the correlation coefficients were the strongest (R = 0.53, p = 0.001, and R = 0.55, p < 0.001, respectively).

Computed Tomography Imaging Shows Airway Wall Thickening in Asthma Patients with Persistent Airflow Limitation
CT imaging revealed a higher RB1 WAR and WTR in subjects with persistent airflow limitation and a tendency towards an increase in those variables in RB10 (Table 3). That group of patients also had lower lumen and airway diameters and decreased lumen areas in RB10 than the remaining individuals (Table 3). Males were characterized by higher airway and lumen areas; however, WAR and WTR were not associated with gender. Additionally, lumen and wall diameters in RB10 remained in a weak positive relationship with height (R = 0.37, p = 0.02, and R = 0.32, p = 0.049, respectively).

Reticular Basement Membrane Thickness Is Not Increased in Asthma Patients with Persistent Airflow Limitation
Surprisingly, RBM thickness was not significantly increased in asthma subjects with persistent airflow limitation ( Figure 2). RBM thickness was also not related to age, sex, asthma duration, and smoking history (data not shown). Likewise, there was no association with spirometry values, asthma severity, or symptom control scores. Among laboratory parameters, thickness of the RBM layer showed a moderate positive correlation with blood neutrophil count (R = 0.35, p = 0.02) and BAL periostin (R = 0.46, p = 0.04), albeit not with systemic biomarkers of type-2 inflammation (blood eosinophil count: R = −0.25, p = 0.09; serum periostin: R = −0.3, p = 0.2). Additionally, patients with a thicker RBM (median value as a cut-off point (6.08 µm)) were characterized by a higher odds ratio (3.95 (95%CI: 1.75-6.2), p < 0.001) of increased BAL periostin concentrations, defined as values above the cut-off point of 0.95 ng/mL. RBM thickness was not related to lung CT airway geometry parameters, except for a weak positive relationship with RB10 wall thickness (R = 0.3, p = 0.04), after adjustment for confounders (β = 0.39 (95%CI: 0.2-0.6)).

Reticular Basement Membrane Thickness is not Increased in Asthma Patients with Persistent Airflow Limitation
Surprisingly, RBM thickness was not significantly increased in asthma subjects with persistent airflow limitation ( Figure 2).

Increased Reticular Basement Membrane Thickness Is Associated with Fibrosis-and Growth-Promoting Pattern of Gene Expression in Bronchial Epithelium
To search for unique airway transcriptome patterns associated with increased RBM thickness, we stratified asthma patients into those with increased vs. decreased RBM thickness, using the median value as a cut-off point (6.08 µm). This comparison revealed 51 genes, for which expression significantly differed between both analyzed subgroups. Among them, 17 genes differed at the significance level lower than 0.001. These genes are described in Table S1 (Supplementary Materials). The thicker RBM was linked with a gene expression profile promoting cell proliferation and growth (CDK20 [8], MIS12 [9], PDS5B [10]), regulating DNA transcription (ZNF594 [11]), apoptosis (RNF157 [12], HTRA2 [13]), and cargo transporting (SLC39A13 [14], SLC16A6 [15], IPO13 [16]), as well as related to the ECM components (LAMA3 [17], USH2A [18]). Moreover, we documented higher mRNA levels of CLEC9A (C-type lectin domain containing 9A) [19], which encodes an endocytic receptor specialized in the uptake and processing of dead cell material, and OAS2, encoding 2 ,5 -oligoadenylate synthase, engaged in the cell antiviral response [20]. Interestingly, we also identified differentially expressed genes involved in neuronal cell proliferation, growth, and function. For example, increased RBM thickness was associated with the lower expression of BRINP2, which inhibits neuronal cell proliferation (Table S1) [21].
Further correlation analyses identified several genes in which expression was positively (101 genes) or negatively (39 genes) associated with RBM thickness. Among them, 19 and nine, respectively, had particularly strong correlation coefficients (i.e., ≥ 0.6 or ≤ −0.6). This set of genes is listed in Table 4. Table 4. Epithelial cell transcriptome-genes for which expression remained in strong positive or negative associations with reticular basement membrane thickness (correlation coefficients ≥ 0.6 or ≤ −0.6).

MIS12
MIS12 kinetochore complex component 0.745 Plays an essential role in chromosome segregation [9]. The higher expression suggests increased proliferation potency of airway epithelial cells.

NEK5
NIMA related kinase 5 0.625 Regulates centrosome integrity and is involved in cell proliferation [24]. The higher expression suggests increased proliferation potency of airway epithelial cells.

ORC5
The origin recognition complex subunit 5 0.603 Promotes DNA replication [25]. The higher expression might indicate increased proliferation potency of airway epithelial cells.

AMBRA1
Autophagy and beclin 1 regulator 1 0.626 Regulates autophagy, cell survival, and proliferation [26]. A positive correlation with reticular basement membrane (RBM) thickness suggests that autophagy mechanisms are likely involved in airway remodeling pathology.

GAD1
Glutamate decarboxylase 1 0.696 An enzyme catalyzing gamma-aminobutyric acid (GABA) production. GABA may be produced by human epithelial cells, likely relaxing airway smooth muscle [27]. Higher expression of GAD1 in airway transcriptome has also been documented in chronic obstructive pulmonary disease (COPD) patients [28].

MYO1B
Myosin I B 0.658 A motor protein related to the actin filament-based cell organization and movement [29]. Higher airway epithelial cell expression might reflect epithelial-to-mesenchymal transition, a structural modification, whereby epithelial cells lose cell-cell polarity and display mesenchymal features, such as migratory properties [5].
An actin-binding protein involved in cell cycle progression, signal transduction, cell motility, apoptosis, and gene expression [30]. Higher expression in airway epithelium may also reflect the epithelial-to-mesenchymal transition.

TPM1
Tropomyosin alpha-1 chain 0.603 An actin-binding protein involved in the contractile cell system [31]. The higher expression may also reflect the epithelial-to-mesenchymal transition.

Otoferlin 0.612
A protein involved in the vesicular release. Pathological mutations in OTOF are related to deafness. Although its exact role in asthma pathology is unknown, it may be involved in airway mucus secretion [32].

TPH1
Tryptophan hydrolase 1 0.643 A hydroxylase that catalyzes the first rate-limiting step in the biosynthesis of serotonin [33], an important inflammatory mediator in asthma. Serotonin modulates the function of various inflammatory cells, such as mast cells, eosinophils, dendritic cells, and lung epithelium.
Higher levels have been shown in asthma patients' airways; however, its exact source in the lungs is unknown [34].

RNF34
Ring finger protein 34 0.683 A ring finger containing protein, involved in protein-protein and protein-DNA interactions. It likely protects epithelial cells from premature apoptosis [35].

EEF2
Eukaryotic translation elongation factor 2 0.658 An essential factor in cell protein biosynthesis. The higher expression suggests increased protein translation potency in asthma patients with thicker RBM [36].

STAP1
Signal transducing adaptor family member 1 0.609 Positively regulates gene expression and signal transduction. The exact role in asthma pathology is unknown; however, it is abundant in activated human B cells [37].
An interferon-induced antiviral ribonuclease L that destabilizes double-stranded viral RNA. It plays a critical role in the cellular innate antiviral response. A positive correlation with RBM thickness suggests that viral infections, essential asthma exacerbation triggers, might also be involved in airway remodeling pathology. The encoded protein also regulates cell apoptosis, growth, differentiation, and gene expression [20].

LAMA3
Laminin subunit α 3 0.615 A protein belonging to the laminin family. Laminins are essential for the formation and function of the basement membrane [17]. Higher airway expression has been linked with COPD and idiopathic pulmonary fibrosis patients [38].

PCDHA3
Protocadherin alpha 3 0.609 A neural cadherin-like cell adhesion protein. Its role in airways is unknown; however, it might be related to the neural network and synaptic organization in airway remodeling pathology [39].

PDS5
Cohesin associated factor B −0.643 A protein that negatively regulates cell proliferation [10]. The lower expression suggests increased proliferation potency of airway epithelial cells.

RNF13
Ring finger protein 13 −0.636 A ring zinc finger containing protein, promoting cell apoptosis [43]. The lower expression might indicate the anti-apoptotic potency of airway epithelial cells.

MRPL37
Mitochondrial ribosomal protein L37 −0.66 A large subunit of mitochondrial ribosomal proteins, which catalyzes protein production within the mitochondrion. Recent investigations have uncovered this protein's pro-apoptotic role [44,45]; therefore, its lower expression might indicate anti-apoptotic potency of airway epithelial cells.

CACNA1G
Calcium voltage-gated channel subunit α1 G −0.615 A transmembrane voltage-dependent calcium channel promoting cell death. The decreased expression suggests lower apoptotic activity [46]. A stretch-gated ion channel that senses airway stretches [47]. Lower expression in those with a thicker RBM might indicate a weaker response to epithelial cells' mechanical stimuli. Lower PIEZO2 expression has been shown in COPD patients [28].

CA3
Carbonic anhydrase 3 −0.606 A carbonic anhydrase that facilitates the reversible hydrating of carbon dioxide. Lower expression in those with thicker RBM might indicate a possible weaker contraction ability or a worse response to oxidative stress in airway epithelium because CA3 is upregulated by HIF1-α [48].

KIF11
Kinesin family member 11 −0.645 A motor cell protein involved i.a. in cilia beating and cell division. The lower expression might suggest impaired cilia function in those with a thicker RBM [49].

C21orf58
Chromosome 21 open reading frame 58 −0.629 Exact biological roles are unknown In Figure 3, we present the cross-correlation matrix and hierarchical clustering (a heatmap) of 28 genes, for which expression correlated best with RBM thickness, i.e., with correlation coefficients ≥ 0.6 or ≤ −0.6. This Figure demonstrates which genes (among those 28) were associated well with each other. For example, some of the genes for which expression remained in a strong inverse relationship with RBM thickness, e.g., MRPL37, were also negatively associated with the mRNA level of genes promoting cell proliferation and growth, but determining RBM thickness positively (e.g., ORC5, CDK20, EEF2). On the other hand, a strong association was demonstrated between mRNA levels of TPMI1 and CDK20 or TACC2 and AMBRA1, likely indicating similar regulatory mechanisms of their expressions.
The remaining genes significantly correlated with RBM thickness (i.e., with less strong coefficients) are listed in Supplementary Materials (Tables S2 and S3).
As expected, correlation analysis confirmed most of the genes differentially expressed in patients with higher RBM thickness (e.g., LAMA3 [17], USH2A [18], MIS12 [9], CDK20 [8]). However, we identified other genes potentially contributing to RBM depositions, including TACC2, encoding centrosome-and microtubule-interacting protein [22], NEK5 regulating centrosome integrity [24], ORC5 modulating DNA replication [25], AMBRA1 involved in autophagy, cell survival, and proliferation [26], as well as MYO1B and CORO2A, both controlling various cellular processes, such as cell cycle progression, signal transduction, cell motility, apoptosis, and gene expression [29,30]. Interestingly, many of the genes positively correlated with RBM thickness have anti-apoptotic functions (e.g., RNF34 [35]), are involved in immune responses (OAS2 [20], STAP1 [37]), the cell contractile system (TPM1, TPM2, MYO1B, CORO2A [29][30][31]), neuronal growth and function (PCDHA3 [50]), or have essential enzymatic activities (GAD1 [27], TPH1 [33], EEF2 [36]). On the other hand, several genes inversely associated with the RBM thickness possess opposite biological functions, such as PDS5 [10] which negatively regulates cell proliferation, and apoptosis promoting genes, including RNF13 [43], MRPL37 [44,45], and CACNA1G [46]. Overall, most of the genes associated with RBM thickness could be assigned to ontology groups linked with a cell growth-promoting phenotype (e.g., cell proliferation and metabolism, anti-apoptotic) and fibrosis/remodeling-promoting phenotype (e.g., ECM components and proteins regulating neovascularization) as summarized in Figure 4. This Figure was prepared based on gene ontology terms presented in Tables S4-S6 in Supplementary Materials and the literature data. In contrast, in Figure 5, we present a network of lungspecific interactions of those genes. As has been shown, the processes modulated by the indicated genes have a wide range of interdependencies. We identified several hubs with interlinking genes, associated with basement membrane structure or functions, either directly or through first-degree interactors. By betweenness centrality within the network, HSP90AB1, encoding multifunctional chaperon protein, BIRC3 regulating apoptosis, HPRT1 referred to the purine salvage pathway uric acid production, and TCEA2, involved in RNA polymerase II functioning, formed the largest hubs. The lower nodes referred to the PML gene, modulating DNA transcription, cell apoptosis, and senescence; MYOB1, engaged in actin filament-based movement, and MUC1 or CLEC7A encoding membrane receptors. The remaining genes significantly correlated with RBM thickness (i.e., with less strong coefficients) are listed in Supplementary Materials (Tables S2 and S3).
As expected, correlation analysis confirmed most of the genes differentially expressed in patients with higher RBM thickness (e.g., LAMA3 [17], USH2A [18], MIS12 [9], CDK20 [8]). However, we identified other genes potentially contributing to RBM depositions, including TACC2, encoding centrosome-and microtubule-interacting protein [22], NEK5 regulating centrosome integrity [24], ORC5 modulating DNA replication [25], AM-BRA1 involved in autophagy, cell survival, and proliferation [26], as well as MYO1B and CORO2A, both controlling various cellular processes, such as cell cycle progression, signal transduction, cell motility, apoptosis, and gene expression [29,30]. Interestingly, many of the genes positively correlated with RBM thickness have anti-apoptotic functions (e.g.,  CACNA1G [46]. Overall, most of the genes associated with RBM thickness could be assigned to ontology groups linked with a cell growth-promoting phenotype (e.g., cell proliferation and metabolism, anti-apoptotic) and fibrosis/remodeling-promoting phenotype (e.g., ECM components and proteins regulating neovascularization) as summarized in  In contrast, in Figure 5, we present a network of lung-specific interactions of those genes. As has been shown, the processes modulated by the indicated genes have a wide range of interdependencies. We identified several hubs with interlinking genes, associated with basement membrane structure or functions, either directly or through first-degree interactors. By betweenness centrality within the network, HSP90AB1, encoding multifunctional chaperon protein, BIRC3 regulating apoptosis, HPRT1 referred to the purine salvage pathway uric acid production, and TCEA2, involved in RNA polymerase II functioning, formed the largest hubs. The lower nodes referred to the PML gene, modulating DNA transcription, cell apoptosis, and senescence; MYOB1, engaged in actin filamentbased movement, and MUC1 or CLEC7A encoding membrane receptors.

Cluster Analysis Integrating Molecular and Remodeling Data Reveals Subgroups of Asthma Patients with Unique Bronchial Epithelial Gene Expression Patterns
A cluster analysis based on epithelial expression levels of genes that correlated best with RBM thickness (28 genes included) and clinical data revealed three different clusters of patients. Clusters one (n = 11) and two (n = 10) were characterized by altered expression of six genes related to cell growth and proliferation or possessing anti-apoptotic properties (i.e., increased expression of CDK20, RNF34, GAD1, PDS5 and lower of RNF13 and MRPL37), involved in ECM composition (LAMA3), innate immune response (OAS2), contractile system (TPMI1, PIEZO2), and TPH1, as compared to cluster three (n = 19). Cluster two had the highest expression of additional genes engaged in ECM structure (USH2A) and cell growth (EEF2, TACC2, NEK5), together with the most significant decline in HPRT, compared with the two remaining collections. Cluster two has also been referred to as having increased expression of the other two genes involved in cell growth (ORC5, AMBRA1) and STAP1 and OTOF, compared to cluster three.
Interestingly, clusters one and two had thicker RBMs in comparison to cluster three (p < 0.001, both), while cluster two referred to the highest values of RB1 and RB10 wall thickness, as compared to the remaining patient collections (p < 0.008 and p = 0.049, respectively). Cluster one had a higher BAL periostin than cluster three (p = 0.01). Demographic and clinical parameters and spirometry values did not differ among the three separated sets of patients.

Cluster Analysis Integrating Molecular and Remodeling Data Reveals Subgroups of Asthma Patients with Unique Bronchial Epithelial Gene Expression Patterns
A cluster analysis based on epithelial expression levels of genes that correlated best with RBM thickness (28 genes included) and clinical data revealed three different clusters of patients. Clusters one (n = 11) and two (n = 10) were characterized by altered expression of six genes related to cell growth and proliferation or possessing anti-apoptotic properties (i.e., increased expression of CDK20, RNF34, GAD1, PDS5 and lower of RNF13 and MRPL37), involved in ECM composition (LAMA3), innate immune response (OAS2), contractile system (TPMI1, PIEZO2), and TPH1, as compared to cluster three (n = 19). Cluster two had the highest expression of additional genes engaged in ECM structure (USH2A) and cell growth (EEF2, TACC2, NEK5), together with the most significant decline in HPRT, compared with the two remaining collections. Cluster two has also been referred to as having increased expression of the other two genes involved in cell growth (ORC5, AM-BRA1) and STAP1 and OTOF, compared to cluster three.

Discussion
In this study, we compared several airway remodeling indices at the large bronchi and peripheral level, including CT airway cross-sectional geometry and histological examination of the mucosa (RBM thickness) in asthma subjects, with a particular focus on individuals with persistent airflow limitation vs. those with normal spirometry after bronchodilator. We also investigated the associations between the bronchial epithelial cell transcriptome and airway structural parameters. We have demonstrated that altered gene expression in the bronchial epithelium may determine RBM thickness, which was surprisingly weakly related to other airway remodeling markers. First of all, the RBM layer was not increased in patients with persistent airflow limitation. It either did not correlate with CT bronchial wall remodeling indices, except for a weak relationship with RB10 wall thickness. On the other hand, we have shown a clear association between the degree of airflow limitation and airway remodeling indices in lung imaging, particularly within the RB1. Therefore, we may speculate that RBM thickening likely progresses independently of the parameters linked with airway geometry and function, such as the rate of the decrease in spirometry values and the altered bronchial wall in CT scans.
Interestingly, airway and systemic biomarkers of type-2 inflammation (e.g., blood and BAL eosinophilia) were higher in patients with persistent airflow limitation and those with thicker walls in CT. In turn, RBM thickness was related to the BAL periostin and blood neutrophilia, albeit not to the systemic type-2 biomarkers. Therefore, we might speculate that separate mechanisms may drive airway remodeling, depending on whether they act at the bronchial mucosa level or airway wall structure. However, type-2 inflammation in the airways is likely involved in both of them.
Another interesting finding of our study is that the RBM thickening did not depend on asthma duration. Therefore, one might speculate that RBM thickening occurs early during the disease, e.g., as a response to the ongoing airway inflammation, and does not progress further due to implemented anti-inflammatory treatment. Noteworthy, this observation mirrors the study performed by Payne et al. [51]. They demonstrated that RBM thickening is already present in young children (median age 13 years) to a similar extent to milder adult asthmatics. These authors did not find an association with patients' age, asthma duration, or lung function, either. Furthermore, it has been demonstrated that even children younger than five years old with confirmed wheezing episodes may present a thicker RBM layer [52].
Epithelial cells are the primary source of BAL periostin, stimulating TGF-β release, which activates the underlying fibroblasts to produce ECM components [6]. Although TGF-β is also secreted by eosinophils [4], the positive relationship between RBM thickness and BAL periostin suggests that the bronchial epithelium may indeed orchestrate mucosal and submucosa changes in airway remodeling pathology. However, the biological role of the thicker RBM in asthma remains uncertain. It has been hypothesized that it constitutes a secondary barrier that compensates for the increased epithelial fragility and leakiness, thus decreasing the entry of harmful particles into the submucosa [6]. However, Rijt et al. [53], based on murine experiments, have shown contrary results, indicating that basement membrane remodeling might promote new sensitizations by causing persistent activation of dendritic cells. Therefore, therapeutic intervention to stabilize RBM thickness, e.g., kinase inhibitors, might be beneficial for severe asthma patients. However, to recommend such a therapeutic approach, the specificity of RBM thickening and its determinants need to be comprehensively investigated.
In our study, analysis of epithelial cell transcriptome revealed several gene groups associated with RBM thickness. As expected, in patients with thicker RBM, we demonstrated a higher expression of genes related to ECM proteins, such as LAMA3 or USH2A. They encode the laminin subunit α 3 and usherin, respectively, proteins that are essential for the basement membrane structure and function, widely expressed in many organs, including the lungs [17,18]. Several epithelial-mesenchymal growth regulators may increase LAMA3 expression, such as epidermal growth factor (EGF) [17]. In turn, usherin contains laminin EGF-like and many fibronectin motifs, and thus, might be upregulated by laminin [18]. Therefore, both proteins are likely involved in increased ECM deposition and could participate in airway remodeling pathology. In this context, the decrease in expression of LTBP1 can also be essential. This gene encodes a protein that maintains TGF-β in the latent ECM bound form [54]. Therefore, its downregulation might contribute to the increased TGF-β activity, promoting ECM remodeling and fibrosis.
Interestingly, RBM thickness remained in a strong positive association with the higher expression of multiple genes encoding proteins involved in cell activation, proliferation, and growth or with an anti-apoptotic function. On the other hand, some genes stronglybut inversely-related to RBM thickness have opposite roles, i.e., either antiproliferative or pro-apoptotic properties. This observation suggests that the airway epithelial transcriptome linked with thicker RBM promotes cell proliferation and growth, further reinforced by the upregulated genes referring to the increased cell metabolism, protein synthesis, cargo transport, and nucleic acid processing.
Another important finding of our study is the lower expression of genes involved in the formation of tight cell junctions. This outcome is in line with other reports pointing to the disrupted airway epithelial barrier in asthma [7]. It is believed that cleaving tight junctions results from external factors, such as allergens, pathogens, air pollutants, and an intrinsic predisposition [7]. Our report suggests that increased permeability of the epithelium may also contribute to subepithelial changes and airway remodeling. However, it is worth noting that parallel to the decreased cadherin mRNA levels, we also documented higher expression of PNN, an E-cadherin transcription activator, likely serving as a compensatory mechanism [55].
An airway epithelial cell transcriptome profile associated with thicker RBM might also reflect epithelial-to-mesenchymal transition, a process whereby epithelial cells acquire characteristics of fibroblasts or mucosa cells [6]. During such a structural modification, epithelial cells lose cell-cell polarity and adhesion, to become migratory and develop mesenchymal features [5]. For example, this is reflected in our data by the increased expression of actin-related proteins, such as TPM1, TPM2, MYO1B, and CORO2A.
The next group of genes associated with a thicker RBM, which merits comment, contains those referred to as the oxidative stress response genes. Some of these genes were reported as upregulated, such as HSP90AB1, promoting cell adaptation to environmental changes [56]; SCARA3, the product of which depletes reactive oxygen species [57], and DNMT3A, encoding an enzyme involved in CpG or non-CpG motifs DNA methylation [58]. However, decreased expression of CA3 [48] and HPRT [40], both induced by hypoxiainducible factor 1α (HIF1-α), a transcription factor that adapts the cell during oxygen deprivation, might suggest an impaired adaptation to oxidative stress in airway remodeling pathology. Furthermore, downregulated HPRT might be related to the local overproduction of uric acid [40], likely stimulating cell proliferation and new blood vessel formation [41,42]. Increased neovascularization might also be related to upregulated FGD5 [59]. Interestingly, bronchial epithelial cells from patients with a thicker RBM also showed a gene profile that promotes nerve regeneration and functioning; therefore, it seems that neural network development may also occur in asthma remodeling pathology.
Another intriguing association shown in our data refers to the significant positive relationship between RBM thickness and the mRNA level of the gene encoding a decarboxylase, which catalyzes γ-aminobutyric acid (GABA) production. It has been recently demonstrated that the airway epithelium is a predominant source of endogenous GABA, contributing to airway smooth muscle tone relaxation [27]. Higher GAD1 expression in those with thicker RBM likely indicates susceptibility to increased GABA epithelial production, partially explaining the lack of association between RBM thickness and spirometry values. Interestingly, the higher expression of this gene in the airway transcriptome has also been documented in chronic obstructive pulmonary disease (COPD) patients [28].
The last important issue that merits comment, is the relationship between the RBM layer and the expression of genes contributing to the immune response. A thicker RBM is strongly determined by upregulated STAP1, encoding a protein abundant in the stimulated B cells [37] and, to a lesser extent, TRBV11-1, encoding the variable domain of the T cell receptor [60]. Interestingly, a thicker RBM was also positively associated with gene transcripts encoding the proteins involved in the antiviral response and other innate immunity elements, such as AMBRA1, involved in autophagy regulation [26], and MUC1 protecting cells against bacterial, viral, or enzyme attack [61].
Our study has several limitations. First of all, we did not analyze the airway transcriptome in healthy controls for ethical reasons. Therefore, we cannot exclude that similar genes also might determine the RBM thickness in controls. Asthma patients enrolled in this study were relatively old and had mostly moderate to severe disease; thus, the findings may not represent younger subjects with milder disease. We investigated bronchial brush biopsy; therefore, we speculate that the outcomes mainly refer to the airway epithelium and to a lesser extent, lung infiltrating inflammatory cells. We analyzed airway transcription profiles at a single time point; therefore, we cannot exclude their changes in time. Statistical associations reported here may not necessarily indicate cause-effect relationships. Finally, the clinical relevance of the demonstrated associations in the airway transcriptional profile in terms of disease severity, progression, and relation to airway remodeling requires further investigation.

Patients
Our patient group consists of 40 non-smoking moderate to severe asthma subjects aged 30-65 years, with at least a 3-year history of confirmed asthma and no exacerbation within six months before the enrollment. About half of the patients (n = 21) were characterized by persistent airflow limitation; the remaining had normal spirometry before or after the bronchodilator. Except for biological treatment, all asthma medications were permitted, including oral corticosteroids at a daily dose equivalent to ≤ 10 mg of prednisolone, unless the amount was the same during the preceding three months.
Diagnosis of asthma was established by a physician, based on recurrent respiratory symptoms (wheeze, cough, shortness of breath, and chest tightness), together with current or historically documented post-bronchodilator increase in FEV 1 of at least 200 mL and 12% from the baseline [62]. Persistent airflow limitation was defined as a FEV 1 /VC index below 0.7 or FEV1 lower than 0.8 of predicted value after bronchodilator. Spirometry and bronchial reversibility test were performed according to the American Thoracic Society standards [63], using a Jaeger MasterLab spirometer (Jaeger-Toennies GmbH; Hochberg, Germany).
Ex-smokers were eligible for the study if they stopped smoking at least five years ago, with a history of fewer than seven pack-years. Arterial hypertension, diabetes mellitus, and hypercholesterolemia were allowed in study participants.
Definitions of asthma severity and co-morbidities, as well as other exclusion criteria, have been provided in Supplementary Materials.
The study was conducted following the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of the Jagiellonian University (approval number: KBET/151/B/2013). All subjects gave written informed consent to participate in the study.

Lung CT
Lung CT with an assessment of airway cross-sectional geometry was performed one hour after 400 µg albuterol administration using 64-raw multidetector computed tomography (Aquilion TSX-101A, Toshiba Medical Systems Corporation, Otawara, Japan) in a helical scanning mode, without administration of intravenous contrast medium. Automated program AW Server (General Electric Healthcare, Wauwatosa, WI, USA) was used to quantify the airway cross-sectional geometry in the RB1 and RB10; including an average diameter of the lumen and airway, average wall thickness, lumen, and wall area, WAR and WTR, both defined in Supplementary Materials.

Bronchofiberoscopy, Endobronchial Biopsy, and Bronchial Brush Biopsy
Bronchofiberoscopy was performed according to the guidelines of the American Thoracic Society [64] using the bronchofiberoscope BF 1T180 (Olympus, Tokyo, Japan) with local anesthesia (2% lidocaine) and in mild sedation (0.05-0.1 mg fentanyl and 2.5-5 mg midazolam given intravenously). During that procedure, 2-3 endobronchial biopsies were taken from the right lower lobe (the carina between B9 and B10) together with the brush biopsy. Collected endobronchial specimens were immediately fixed in 10% neutral buffered formalin solution (Sigma-Aldrich, Saint Luis, MO, USA) and sent to the Pathology Department for further analysis. Brushes were immediately immersed in TRIzol Reagent (Thermo Fisher Scientific, Carlsbad, CA, USA).

Bronchoalveolar Lavage Fluid Analysis
The cytospin preparations (Thermo Scientific, Walthman, MA, USA) were made from bronchoalveolar lavage fluid (BALF) samples and stained with the May-Grunwald Giemsa dye. A total of 1000 cells in each preparation were counted; the results were shown as a percentage of all inflammatory cells (with the exception of epithelial cells). The remaining BALF sample was centrifuged; the supernatant was frozen in aliquots and stored at −70 • C until analysis.

Blood Laboratory Investigations
Fasting blood samples were drawn from the antecubital vein using minimal stasis between 8:00 and 11:00 A.M. Complete blood cell and platelet count and fibrinogen were assayed by routine laboratory techniques. High-sensitivity C-reactive protein and IgE were measured by latex nephelometry (Siemens, Marburg, Germany). Blood samples were drawn into serum tubes and centrifuged, similarly to BALF at 2000× g for 20 min. The supernatant was frozen in aliquots and stored at −70 • C until analysis.

Histologic Examination
Tissue specimens were processed routinely, as described in our previous publication [65]-briefly, the 2 µm paraffin-embedded sections were cut and stained by hematoxylin and eosin. The slides were photographed by a Nikon D5300 camera attached to the Zeiss Axioscope microscope with a 100× oil immersion lens. The images were analyzed by the AnalySIS 3.2 software (Soft Imaging System GmbH, Germany). The RBM thickness was measured along the biopsy's epithelial surface (Figure 2) according to the orthogonal intercept method suggested by Ferrando et al. [66] using arbitrary-distance units. We only analyzed sections covered by non-tangentially cut and well-preserved epithelial cells. For each patient, at least 30 individual RBM measurements were evaluated at intervals of 9.5 µm. The results were expressed as harmonic mean, defined in our previous publication [65].

RNA Preparation for Microarray Analysis
Total RNA was extracted from bronchial brush biopsy samples with TRIzol (Thermo Fisher Scientific, Carlsbad, CA, USA), fractioned by gravity gradient, and isolated using chromatographic columns (A&A Biotechnology, Gdynia, Poland), followed by storage in 96% ethanol at −80 • C. The quality of each RNA sample was assessed using Qiagen QIAxel system (Qiagen, Hilden, Germany), and RNA integrity was verified by agarose gel electrophoresis. RNA was reverse-transcribed into the cDNA library using Syngen Universal Script Reverse Transcriptase (Syngen, Wroclaw, Poland). The product was purified by Syngen PCR ME Mini Kit (Syngen, Wroclaw, Poland) and fluorescently labeled and purified using Kreatech ULS Platinum Bright Red/Orange Kit (Kreatech, Amsterdam, Netherland). Hybridization of cDNA to the Human Genomic 49K Mi ReadyArray (HEEBO) (Microarray Inc., Huntsville, AL, USA) at 37 • C for 24 h was performed, and microarrays subsequently were washed.

Microarray Data Retrieval and Bioconductor
Microarrays were scanned with an InnoScan 900 Microarray Scanner, and hybridization signals were detected using Mapix software (v.6.0.1; Innopsys, Carbonne, France). Since the Innopsys Mapix software failed to perform microarray gridding correctly in some cases, we developed a custom gridding algorithm to refine the results. Moreover, visible artifacts and noise render the task challenging. Therefore, the method based on maximally stable extremal regions (MSER) was applied to coordinate positioning markers and establish the sub-grid lines [67]. Computed results were verified manually and corrected if required. For our study, overall gridding performance exceeded that obtained with the commercial software, ultimately resulting in expression profiles of 33,519 gene products. Microarray data preprocessing (background correction and normalization) was done by Limma R software using the established methods [68,69]. We analyzed genomic data with Bioconductor v.3.7. the software of the R environment v.3.5.0.
Biological processes, molecular functions, and cellular components gene ontology terms associated with genes which expression levels correlated significantly with RBM thickness were derived from GO Biological Process 2018 by Gene Ontology Consortium. Only terms with p-value < 0.05 were shown (Supplementary Materials Tables S4-S6).
Semi-automated analytics platforms Cytoscape 3.7.1 and NetworkAnalyst 3.0 were used to create an interactions network, where individual molecular entities formed the nodes and their interactions the edges.

Statistical Analysis
Statistical analysis was performed with Statistica TIBCO 13.3 and R (version 3.6.1) software.
We used the Shapiro-Wilk test to verify data distribution. Continuous variables, mostly non-normally distributed, were reported as a median with interquartile range or as a mean with standard deviation, as appropriate. They were compared by the Mann-Whitney U-test or unpaired t-test, respectively. Categorical variables were given as percentages and compared by χ2 test. Lung CT variables and RBM thickness, were Box-Cox transformed and a one-way covariance analysis (ANCOVA) was performed to adjust for potential confounders, including age, sex, and BMI. To evaluate the relationship between continuous variables, a Spearman rank correlation test or univariate linear regression model with adjustment for age, sex, and BMI (CT parameters and RBM thickness) were performed. To calculate odds ratio with 95% confidence interval (CI), the cut-off point of BAL periostin concentration was calculated based on the receiver operating characteristic (ROC) curve.
We analyzed the Pearson linear and the Spearman rank correlation coefficients to study relationships between expressions of all genes and RBM thickness. The significance levels were established as p < 0.001 and p < 0.05, for the Pearson and Spearman rank correlation coefficients, respectively. We pointed out two genes group. The first was composed of 28 genes characterized by a powerful impact on the RBM thickness (at least 36%) with the absolute value of the correlation coefficient 0.6 and more, or −0.6 and less. The second group consists of genes with the correlation coefficients 0.4 to 0.599 and −0.599 to −0.4. All data regarding gene analyses, such as correlation tests and gene expression comparisons (Mann-Whitney U-test), were adjusted using the Benjamini-Hochberg correction.
A cluster analysis was performed using the k-means clustering method, based on mRNA levels of genes correlated best with RBM thickness. We obtained three different clusters compared by the covariance analysis (ANCOVA), Kruskal-Wallis-test, or χ2 test, as appropriate.
Results that presented a p-value less than 0.05 were considered statistically significant unless otherwise stated in the text.

Conclusions
Our data showed that RBM thickening likely progresses independently of the parameters linked with airway geometry and functioning, such as the rate of the decrease in spirometry values and the thickening of the bronchial walls in CT scans. A thicker RBM was associated with the airway epithelial cell transcriptome promoting cell activation, proliferation, and growth. Increased RBM thickness was also linked with altered expression of genes likely involved in the structural remodeling of the airways, such as increased ECM depositions, promoting neovascularization, and nerve functioning. However, the clinical relevance of the presented relationships, in terms of asthma severity, progress, and airway remodeling of larger airways requires further observational and experimental investigations.
Supplementary Materials: This article involves an online Supplement https://www.mdpi.com/ 1422-0067/22/3/998/s1. Table S1: Epithelial cell transcriptome. Asthma patients were divided based on the median value of the reticular basement membrane (RBM) thickness as a cut-off point; genes with higher or lower expression in those with thicker RBM are presented, Table S2: Genes for which expression remained in weaker, albeit significant positive associations with reticular basement membrane thickness (correlation coefficients 0.4 to 0.599), Table S3: Genes for which expression remained in weaker, albeit significant negative associations with reticular basement membrane thickness (correlation coefficients −0.599 to −0.4), Table S4. Biological processes gene ontology terms associated with genes which expression level correlated significantly with reticular basement membrane thickness, Table S5. Molecular functions gene ontology terms associated with genes which expression levels correlated significantly with reticular basement membrane thickness, Table S6. Cellular components gene ontology terms associated with genes which expression levels correlated significantly with reticular basement membrane thickness.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of the Jagiellonian University approved the study (approval number: KBET/151/B/2013).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to patients' origin.

Conflicts of Interest:
The authors declare no conflict of interest.