Comparative study of epithelial entry gene expression of SARS-CoV-2 and other human viral species in asthma: differences by sex, airway location and disease endotype

Epithelial characteristics underlying the differential susceptibility of chronic asthma to SARS-CoV-2 (COVID-19) and other viral infections are currently unclear. By revisiting transcriptomic data from patients with Th2 low versus Th2 high asthma, as well as mild, moderate and severe asthmatics, we characterized the changes in expression of human coronavirus and influenza viral entry genes relative to sex, airway location and disease endotype. We found sexual dimorphism in expression of COVID-19 genes ACE2, TMPRSS2, TMPRSS4, and SLC6A19. ACE2 receptor downregulation occurred specifically in females in Th2 high asthma, while proteases broadly assisting coronavirus and influenza viral entry, TMPRSS2 and TMPRSS4, were highly upregulated in both sexes. Overall, changes in COVID-19 gene expression were specific to Th2 high molecular endotype of asthma, and different by asthma severity and airway location. The downregulation of ACE2 (COVID-19, SARS) and ANPEP (HCoV-229E) viral receptors correlated with loss of club and ciliated cells in Th2 high asthma, while the increase in DPP4 (MERS-CoV), ST3GAL4, and ST6GAL1 (influenza) associated with an increase in goblet and basal activated cells. Overall, this study elucidates sex, airway location, disease endotype and changes in epithelial heterogeneity as factors underlying asthmatic susceptibility, or lack thereof, to COVID-19.


INTRODUCTION
As the COVID-19 pandemic continues to unfold worldwide, new insights are coming to light about the mechanisms and epidemiology of this disease. One of the most puzzling recent observations is an apparent lack of association of SARS-CoV-2 infection with asthma. Reports from Wuhan, China revealed that asthma was not a risk factor for severity and mortality in adult COVID-19 patients -the prevalence of asthma in patients with COVID-19 was only 0.9%, which is substantially lower than that in the adult population of Wuhan (6.4%) 1 . In another study of 140 patients infected with SARS-CoV-2 in Wuhan, none presented with asthma or other allergic diseases 2 . A study by Mahdavinia et al. 3 in a large cohort of adult COVID-19 cases in the United States showed that asthma was significantly associated with longer intubation but was not associated with higher rate of hospitalization, death, or acute respiratory distress syndrome among COVID-19 patients. This is surprising because respiratory viral infections in general are known to associate with asthma exacerbations and more serious illness.
One study proposed an explanation for this by showing that respiratory allergies and controlled allergen exposures were each associated with significant reductions in epithelial expression of ACE2, the entry receptor for SARS-Cov-2 4 . Specifically, ACE2 expression was lowest in individuals with both high levels of allergic sensitization and asthma 4 . Another study confirmed this observation by showing that Type 2 cytokine IL-13 decreases expression of ACE2 in bronchial epithelium, thus providing one mechanistic link to Th2 inflammation 5 . In this study, we pursued the hypothesis that susceptibility of asthmatic epithelium to different viral species depends on differential expression of receptor entry genes in different sexes, airway locales and asthma endotypes.

Transcriptomic data
In this study we used the following two independent gene expression datasets downloaded from Gene Expression Omnibus (GEO, NCBI): (1) GSE67472, providing data for gene expression in bronchial epithelial brushings from healthy controls and subjects with Th2 low and Th2 high asthma and (2) GSE41861, providing epithelial brushing data for mild, moderate and severe asthma in both upper and lower airways. The study comparing Th2 low and Th2 high asthma by Woodruff et al. 6 is the most highly cited study providing basis for molecular endotyping of asthma patients based on patterns of Type 2 inflammation. Clinical descriptions of asthma patients and airway sampling procedures can be found in GEO Datasets and in the original publication of this study 6 . Gene annotations were updated using technology-corresponding annotation libraries in R. Each individual GEO dataset underwent consistent handling, implementing standard contrast fits on Series Matrix File data uploaded to GEO by the original study groups, preserving original normalization protocols to maintain consistency with gene signatures published by original studies.

Clustering of epithelial and viral entry genes
Single cell analysis data used to generate epithelial subset gene list for unsupervised clustering of epithelial and viral genes is from Nature Medicine study by Vieira Braga et al. 7 providing a cellular census of human lung transcriptome in health and asthma. Exact genes and epithelial subset clustering structure in health and asthma can be found in Fig. 3c in the original Vieira Braga's publication. Viral genes were selected based on published works of human coronavirus and influenza entry genes in epithelial cells. Epithelial-viral gene expression data set was then assembled from GSE67472 data and included both males and females, and healthy, Th2 low and Th2 high individuals. The resulting data set was subjected to hierarchical clustering using One Minus Pearson correlation metric and average linkage method. Similarity matrix was derived from clustered data using Pearson correlation metric. All procedures were performed in MORPHEUS, https://software.broadinstitute.org/morpheus (Broad Institute), matrix visualization and analysis software. Clustering structure in Fig. 1d in this publication is nearly an exact match to the above mentioned study, demonstrating alignment in methodology and biological message in regards to epithelial heterogeneity.

Statistical analysis and data visualization
All graphs and statistical analyses were carried out in GraphPad Prism (Graphpad) and Systat 13 (Systat Software, Inc.). Principal Component Analysis (PCA) was executed in PAST 2.17 software using covariance matrices and singular value decomposition. Statistical significance was determined by ANOVA followed by Tukey's post-hoc pairwise testing. All data are represented as mean ± S.E.M.
An alpha level of 0.05 was used as a significance cut-off. Marginal values (0.05>p<0.1) are also shown.

Data availability
All raw data used in this study, including series matrix file data, can be accessed in GEO under the accession numbers GSE67472 and GSE41861. scRNA-seq data that support the findings of this study can be found in the above mentioned study by Vieira Braga et al 7 .

RESULTS AND DISCUSSION
We revisited publicly available transcriptomic data from two independent studies of bronchial and nasal epithelium for an unbiased assessment of epithelial expression of COVID-19, influenza, and other human coronavirus related genes in (1) Th2 low versus Th2 high asthma 6 in bronchial epithelial brushings (GSE67472) and (2) expression of these genes in mild, moderate, and severe asthma in both upper and lower airways (GSE41861). We found that dysregulation of these genes in asthma is both sexually dimorphic and airway location dependent. For SARS-CoV-2 genes, we found marginally significant downregulation of viral receptor ACE2 in females in Th2 high asthma and significant downregulation in males in mild asthma (upper airway) (Fig. 1a). Expression of specific ACE2partnering membrane protein SLC6A19 (B 0 AT1) was also found to be downregulated in males in mild asthma (lower airway) and moderate asthma (upper airway) (Fig. 1a). Interestingly, we found an opposite pattern of expression for two related membrane-bound serine proteases, TMPRSS2 and TMPRSS4, which promote SARS-CoV-2 virus epithelial entry by processing the spike glycoprotein resulting in the release of viral contents into the host cell cytosol. Expression of both proteases was significantly upregulated specifically in Th2 high asthma (Fig. 1a). We also found upregulation of TMPRSS4 in females in severe asthma (lower and upper airway) and TMPRSS2 in males in mild and moderate asthma (lower airway) (Fig. 1a). Expression of another protease suspected in assisting viral entry, FURIN, was not found to be significantly different in any dataset (data not shown). Such differences due to sex deserve further investigation, as COVID-19 resulted in higher mortality and severity in males 8 . Downregulation of ACE2 expression predominantly in mild asthma is also interesting, as asthma is a very heterogeneous disease. The Center for Disease Control (CDC) currently lists moderate to severe asthma as a risk factor for severe COVID-19 illness (last updated on January 7). Whether the downregulation of ACE2 confers a protective effect in mild asthma as opposed to more severe forms of disease remains to be investigated. Another study showed that the use of inhaled corticosteroids in asthma was associated with lower expression of ACE2 and TMPRSS2 9 . Of note, severe asthma patients are the ones most likely to use corticosteroids, and yet we found ACE2 downregulation only in mild asthma. Principal Component Analysis (PCA) analysis of covariance in several human virus entry genes (ANPEP, DPP4, CTSL, CTSB, ST3GAL4, ST6GAL1) and associated viral immunity genes (IRAK3, IDO1, NOS2, OAS1, MX1, TNFSF10) demonstrated that dysregulation of viral immunity in general is characteristic of the Th2 high but not Th2 low subset of asthma patients (Fig. 1b).
We further performed a comparative analysis of SARS-CoV-2 genes with genes promoting the entry of other human coronaviruses and influenza. We found that ANPEP, the entry receptor for HCoV-229E, was significantly downregulated specifically in Th2 high asthma, as well as downregulated in females in mild, moderate, and severe asthma (lower airway) (Fig. 2a).
Interestingly, HCoV-229E appears to be the human coronavirus least associated with asthma exacerbations. Studies have failed to find HCoV-229E in either adult or childhood episodes of asthma exacerbation 10 . However, asthma has been shown to be a risk for MERS-CoV (Middle Eastern Respiratory Syndrome) coronavirus infection. In one study, 31% of those with MERS had asthma 11 . DPP4, the entry receptor for MERS-CoV, was downregulated in females in Th2 low asthma (trending up in Th2 high asthma) and significantly upregulated in females in mild and moderate asthma (lower airway) (Fig. 2a). Two genes that promote the synthesis of receptors for human influenza strains, ST3GAL4 and ST6GAL1, were specifically upregulated in Th2 high inflammation in both sexes (Fig.   2a). SARS-Cov-1 (SARS) also relies on the ACE2 receptor as an entry mechanism. The SARS (SARS-CoV-1) outbreak also did not affect asthma exacerbations in children 12 .  17 . In a single cell study in Macaca mulatta and humans, ACE2 was also detected primarily in ATII cells but also in club cells 18 . Together, this shows that viral gene expression is restricted to specific epithelial cell subsets.
Unfortunately, there are no studies comparing viral genes in regards to epithelial heterogeneity across Th2 low and Th2 high asthma phenotypes. By comparing known markers of epithelial cell subsets between healthy, Th2 low, and Th2 high asthma, we found evidence for a remarkable reshaping of epithelial heterogeneity in Th2 high asthma (Fig. 2b). We found an increase in markers of basal cycling, basal activated, goblet, mucous ciliated cells, and ionocytes in allergic inflammation.
At the same time, there was a significant decrease in markers of ciliated cells, submucosal glands, club (Clara) and tuft cells. This pattern of epithelial dysregulation was highly specific to Th2 high inflammation, confirming at high resolution loss of epithelial differentiation and specialized cell subsets expressing COVID-19 receptor (club cells, ciliated cells). To gain further insights whether expression of viral genes may correspond to epithelial heterogeneity pattern, we clustered top markers describing each epithelial subset (guided by a recently published single cell study characterizing novel airway epithelial cell subsets in asthma 7 ) with entry genes of different viral species in "Th2 high vs. healthy control" gene expression analysis (Fig. 2c). Our analysis suggested that genes corresponding to different viral species may co-express with markers of different epithelial subsets. In particular, ACE2 and ANPEP were co-expressed with markers of club, tuft, ciliated cells, and submucosal glands, all downregulated during Th2 high (but not Th2 low) inflammation (Fig. 2c).
TMPRSS2 expression very closely correlated with mucous ciliated cells, while TMPRSS4, DPP4 and influenza genes correlated with goblet and basal activated cells, all increasing in Th2 high inflammation (Fig. 2c). Although we acknowledge that this is only an exploratory correlation analysis, which requires further rigorous validation at the single cell analysis level, it does show importance of relating viral gene expression to epithelial heterogeneity. Such localization is also supported by single cell studies. This also explains why there is a complete discrepancy in the expression of ACE2 and TMPRSS2/TMPRSS4 in Th2 high inflammation, since these genes differentially follow epithelial subset divergence during inflammation. It is also very important to acknowledge that while ACE2 is a specific to certain human coronavirus species (SARS-CoV-1, SARS-CoV-2, and HCoV-NL63), TMPRSS2 and TMPRSS4 are certainly not. These proteases are important for cleavage of the spike protein of multiple pathogenic viral species, including MERS-CoV 19 as well as influenza strains 20 .
In summary, changes in epithelial specialization during Th2 high inflammation, as well sex, location and disease endotype may underlie differential susceptibility to SARS-CoV-2 and other