Transcriptomic analysis of COVID‑19 lungs and bronchoalveolar lavage fluid samples reveals predominant B cell activation responses to infection.

The outbreak of the 2019 coronavirus disease (named, COVID‑19), caused by the novel SARS‑CoV‑2 virus, represents a worldwide severe threat to public health. It is of the utmost importance to characterize the immune responses against the SARS‑CoV‑2 and the mechanisms of hyperinflammation, in order to design better therapeutic strategies for COVID‑19. In the present study, a transcriptomic analysis was performed to profile the immune signatures in lung and the bronchoalveolar lavage fluid samples from COVID‑19 patients and controls. Our data concordantly revealed increased humoral responses to infection. The elucidation of the host responses to SARS‑CoV‑2 infection may further improve our understanding of COVID‑19 pathogenesis and suggest better therapeutic strategies.

Introduction

The severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) was first isolated at the end of 2019 in China (1-5) and, as of August 3rd, 2020, almost 18 million infected patients and 686,703 deaths have been reported globally (WHO Situation Report-196). However, the actual number of the infected subject is under-estimated and, indeed, a recent meta-analysis performed on 50,155 patients from 41 studies, showed that the pooled percentage of asymptomatic infection is 15.6% (6).

Even if SARS-CoV-2 shares similarities with the other coronaviruses, the higher diffusion rate and the possibility to induce fatal complications, such as severe pneumonia, acute respiratory distress syndrome (ARDS), thrombosis, septic shock and organ failure, make this virus a major public health threat (7-10). Development of COVID-19 complications seems to be dependent on a dramatic release of proinflammatory factors, such as interleukin (IL)-1β, IL-6, IL-8, tumor necrosis factor-α (TNF-α) and CXC-chemokine ligand 10 (CXCL10) and CC-chemokine ligand 2 (CCL2) in the infected lung tissue and other peripheral organs (2,11-13), which ultimately leads to a reaction known as cytokine release syndrome (CRS). It is likely that CRS promotes a self-sustaining inflammatory process that contributes to the respiratory failure and the systemic manifestations observed in COVID-19 patients (14). A multicenter study of 150 confirmed COVID-19 cases in Wuhan, China, identified as predictors of mortality both elevated ferritin (15,16) and IL-6 levels, which strengthen the hypothesis that fatality events may be due to a virus-driven hyperinflammation (2,11).

The rapid worldwide diffusion of SARS-COV-2 has propelled both basic science and clinical research studies for the elucidation of the pathogenetic mechanisms underlying COVID-19. The emerging observation that a significant percentage proportions of individuals are asymptomatic, not only suggests that SARS-CoV-2 may have a longer incubation period and higher transmission rate, as compared to other coronaviruses, but also advocates potential differences in the host immune responses to this virus. It is therefore, of the utmost importance to characterize the immune responses put against SARS-CoV-2 and the mechanisms of hyperinflammation, in order to design better therapeutic strategies for COVID-19. In the present study, we performed a transcriptomic analysis to profile the immune signatures in lung and the bronchoalveolar lavage fluid samples from COVID-19 patients and controls. Our data concordantly revealed increased humoral responses to infection. The elucidation of the host responses to SARS-CoV-2 infection may further improve our understanding of COVID-19 pathogenesis and suggest better therapeutic strategies.

Materials and methods

Dataset selection

The NCBI Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) was interrogated using the terms 'SARS-CoV-2' and 'COVID-19'. The available datasets were shortlisted if: i) they included whole-genome transcriptomic profiling; ii) included human samples; and iii) were not generated on cancer cell lines. Finally, the GSE150316 and the GSE147507 (17) datasets were selected. GSE150316 is a high throughput sequencing dataset of five autopsy samples from patients deceased due to SARS-CoV-2 infection (2-5 technical replicates for each sample were averaged for the downstream analysis) and five negative control samples.

The GSE147507 dataset was generated from three biological replicates of primary human lung epithelium either infected with SARS-CoV-2 (USA-WA1/2020) at a multiplicity of infection (MOI) of 2, for 24 h, or mock infected. Total mRNA libraries were prepared using tTruSeq Stranded mRNA LP and cDNA libraries sequenced using an Illumina NextSeq 500 platform. Raw reads were aligned to the human genome (hg19) using the RNA-Seq Alignment App on Basespace (Illumina). The submitter-supplied pre-processed and normalized gene expression matrix was used for the analysis.

For the transcriptomic analysis of COVID-19 BALF samples, RNA-Seq data from the Genome Sequence Archive of the Beijing Institute of Genomics (BIG) Data Center (https://bigd.big.ac.cn/) (accession no. CRA002390), and from the NCBI SRA database (accession nos. SRR10571724, SRR10571730 and SRR10571732) (18) were used.

Enrichment and network analysis

Functional enrichment analysis was conducted using the web-based utility, Metascape (19). Metascape analysis makes use of public databanks, such as Gene Ontology, KEGG, and MSigDB, and aggregates enriched ontology terms into non-redundant groups, by calculating the similarity between any two terms (19). Metascape uses the hypergeometric test and the Benjamini-Hochberg p-value correction to identify statistically significant enriched terms. Representative terms from the enrichment analysis are presented as a network. Each term is represented by a node, with its size being proportional to the number of input genes belonging to that term, and the color representing its corresponding cluster. Terms with a similarity score >0.3 are linked by an edge. The thicker the edge, the higher the similarity score. The network is visualized using Cytoscape (version 3.1.2) with 'force-directed' layout. One term from each cluster has its description shown as a label.

Computational deconvolution of infiltrating immune cells

In order to evaluate the relative proportions of immune cell subsets in COVID-19 and healthy control samples, we performed a computational deconvolution analysis. To this end, we used the xCell software, a web computational utility that aims at evaluating, by using gene signatures, the relative proportions in a sample of various immune cell types, including immature dendritic cells (iDCs), conventional DCs (cDCs), active DCs (aDCs), plasmacytoid DCs (pDCs), B cells, CD4+ T cells, memory cells, Th1 cells, Th2 and Treg cells and macrophages (20).

Statistical analysis

The differential expression analysis was performed using the DeSeq2 function. The web-based application NeworkAnalist was used for the statistical analyses. Genes with an adjusted P-value <0.05 were identified as differentially expressed genes (DEGs) and selected for further analysis.

Linear regression and Spearman's correlation were performed to compare the expression levels of genes in COVID-19 samples as compared to healthy control samples, the GSE150316 and the GSE147507 datasets.

For the analysis of the deconvolution data, normality was first assessed using the Shapiro-Wilk, D'Agostino-Pearson and Kolmogorov-Smirnov tests. Based on the results, differential analysis as performed using the non-parametric Mann-Whitney U test.

The GraphPad Prism (version 8) software (GraphPad Software, Inc.) and the SPSS software (SPSS, Inc.) were used for the statistical analysis and the generation of the graphs. Unless otherwise stated, P<0.05 was considered to indicate a statistically significant difference.

Results

Network and enrichment analysis of SARS-CoV-2 infection

In order to determine the transcriptomic signature of lung tissues from COVID-19 patients, we analyzed the GSE150316 RNA-Seq dataset. A total of 55 differentially expressed genes was found, of which 32 were upregulated and 23 downregulated. Gene term enrichment analysis identified GO:0002377: immunoglobulin production, GO:0006959: humoral immune response and GO:0002758: innate immune response-activating signal transduction, as significantly enriched among the upregulated genes (Fig. 1A). A heatmap of the genes belonging to the GO:0002377 (immunoglobulin production) category is presented in Fig. 1B. Among the downregulated genes, GO:0071236: cellular response to antibiotic, GO:0048511: rhythmic process, GO:0042698: ovulation cycle, GO:0019221: cytokine-mediated signaling pathway and GO:0001503: ossification, were found to be significantly enriched (Fig. 1A).

Figure 1

Enrichment and network analysis of COVID-19 lung samples. (A) Hierarchical clustering of the top most enriched terms among the genes signifi-cantly upregulated and downregulated in the GSE150316 dataset. The heatmap is colored based on the p-values, and grey cells indicate the lack of significant enrichment. (B) Gene expression heatmap of the genes belonging to the GO:0002377 (immunoglobulin production) category, in lung biopsies of COVID-19 patients and control samples, as determined in the GSE150316 dataset.

We have previously interrogated the GSE147507 dataset, which included transcriptomic data from primary human bronchial epithelial cells infected in vitro with the SARS-CoV-2 virus (18). Here, we compared this gene signature to the transcriptomic signature of lung biopsies from COVID-19 patients. The publicly available GSE150316 dataset was used in order to perform a correlation analysis on the modulation of the genes perturbed upon SARS-CoV-2 infection and the corresponding genes in GSE147507. A total of 9602 genes were in common between the two datasets. As shown in Fig. 1A, a moderate but significant correlation is found in the transcriptomic profile of in vitro infected bronchial epithelial cells and lung samples from patients (Fig. 2A).

Figure 2

Enrichment and network analysis of COVID-19 samples. (A) Scatter plot showing the correlation of the gene expression profile between the GSE150316 and the GSE147507 datasets. (B) Hierarchical clustering of the top most enriched terms among the genes significantly upregulated and downregu-lated in the GSE150316 and the GSE147507 datasets. The heatmap is colored based on the p-values, and grey cells indicate the lack of significant enrichment. (C) Representative terms from the enrichment analysis are presented as a network, visualized as a 'force-directed' layout. Description of each term is shown as a label.

Gene term enrichment analysis for the significantly modulated genes identified pathways in common between the GSE147507 and the GSE150316 datasets. The common enriched terms for the upregulated genes in the two datasets were: 'humoral immune response' (GO:0006959) and 'leukocyte migration' (GO:0050900) (Fig. 2B).

Interestingly, the top terms enriched among the downregulated genes in the GSE150316 dataset were: 'TNF signaling pathway' (hsa04668), 'cytokine-mediated signaling pathway' (GO:0019221), 'myeloid leukocyte activation' (GO:0002274) and 'regulation of cytokine production' (GO:0001817) (Fig. 2B). Representative terms from the enrichment analysis and their functional connections are presented as a network (Fig. 2C).

Deconvolution analysis of infiltrating immune cells in lung samples from COVID-19 patients

We next characterized the relative proportions of infiltrating immune cells in the lungs of COVID-19 patients. A shown in Fig. 3A, a moderate, but not significant, increase in the immune score and microenviroment score was detected for the COVID-19 lung samples. Also, a moderate, non-significant increase in the percentage of infiltrating basophils and aDCs was observed (Fig. 3B). Analysis of the lymphoid cells in the lungs of COVID-19 patients revealed a significant higher proportion of infiltrating B cells upon SARS-CoV-2 infection, along with a moderate, non-significant increase in NKT and Th1 cells (Fig. 3C).

Figure 3

Deconvolution analysis of infiltrating immune cells in lung samples from COVID-19 patients and controls as determined in the GSE150316. (A) Histogram plot for the ImuneScore, StromaScore an MicroenviromentScore in lung samples from COVID-19 patients and controls, as determined in the GSE150316. (B) Relative proportions of infiltrating myeloid cells in lung samples from COVID-19 patients and controls, as determined in the GSE150316. (C) Relative proportions of infiltrating lymphoid cells in lung samples from COVID-19 patients and controls, as determined in the GSE150316.

Characterization of the transcriptomic profile of BALF samples

Next, we compared the gene signature of BALF samples from COVID-19 patients and controls. A total of 3003 genes were found to be modulated in SARS-CoV-2 patients (adjusted P-value <0.05 and ǀfold-changeǀ >2), with 1745 genes being upregulated and 1258 genes downregulated. As shown in Fig. 4A, among both the upregulated and downregulated genes, pathways related to cell morphology (GO:0051017; GO:0030155; GO:0030036; GO:0030031) and survival (GO:0010942; GO:0097190) were significantly enriched.

Figure 4

Characterization of the transcriptomic profile of BALF samples from COVID-19 patients and controls. (A) Hierarchical clustering of the top most enriched terms by genes significantly upregulated and downregulated in the BALF samples of COVID-19 patients as compared to control samples. (B) Hierarchical clustering of the top predicted transcription factors involved in the expression of the genes significantly upregulated and downregulated in the BALF samples of COVID-19 patients as compared to control samples. (C) Deconvolution analysis of infiltrating immune cells in BALF samples from COVID-19 patients and controls.

Analysis of the transcription factors identified RELA, NFKB1, USF2 and SP1, as putatively involved in the regulation of the differentially expressed genes (Fig. 4B).

Immune cell deconvolution analysis revealed a trend of higher proportion in B cells (both naïve, memory and plasma cells), along with an increase in CD4 memory T cells, CD8 T cells and DCs (cDCs, iDCs and pDCs) (Fig. 4C).

Discussion

The characterization of the exact pathogenetic mechanisms by which SARS-CoV-2 induces multiple organ damage are of immediate importance. Emerging data seem to indicate that beside lungs, other organs, including heart, kidney and the central nervous system may also be affected in COVID-19 (21,22). Patients may show proteinuria, hematuria and increased creatinine levels (21), and may suffer from neurological symptoms, such as headache, epilepsy, disturbed consciousness, anosmia and dysgeusia (22). Some COVID-19 patients also develop thromboembolic events, with elevation of D-dimer and other procoagulant parameters (23), which may represent a secondary anti-phospholipid syndrome (APS) (24), as well as other autoimmune diseases (25). Indeed, accumulating case reports show that COVID-19 patients tested positive for anti-CL, anti-b2-GPI autoantibodies (26-29), as well as lupus anticoagulant (27,30,31).

The use of gene expression profiling data has been extensively employed for the identification of novel pathogenic pathways and therapeutic targets (32-36) for several disorders including, autoimmune diseases (37-40) and cancer (41,42). A computational analysis was performed in order to characterize the immune response to SARS-CoV-2 infection. To achieve this, we exploited publicly available RNA-seq data, generated from lung biopsies and BALF samples from COVID-19 patients. Our data from lung and BALF samples concordantly show that B cell responses mainly characterize SARS-CoV-2 infection.

It has been already described that SARS-CoV-2 elicits a robust humoral cell response, with virus-specific IgM, IgG and IgA, and neutralizing IgG antibodies following infection. Seroconversion usually occurs in most COVID-19 patients between one to two weeks after overt symptoms, and antibody titers last for weeks, following virus eradication (43). It seems also that protective B cell memory arises following infection, as a recent study of SARS-CoV-2 infection in rhesus macaques found that animals that had resolved the primary infection were resistant to reinfection one month later (44). Also, independent data show that higher virus-specific antibody titers correlated with greater virus neutralization and are inversely correlated with viral load (43). However, higher titers may be associated with more severe clinical cases (45-47), suggesting that the humoral responses may not be sufficient to protect from severe disease. Up to now, there is no evidence that SARS-CoV-2-induced anti-bodies contribute to some of the pathological features observed in COVID-19 patients. However, this possibility should be taken into consideration in light of the above-mentioned data on secondary APS syndrome in some COVID-19 cases. It has been proposed that antibody-dependent enhancement (ADE), may represent, at least one, of the causes of the CRS (48-51). When the virus infects the body, memory B cells are activated while the activation of naive B cells is inhibited. However, both virus-specific antibodies and antibodies cross-reacting with other similar virus strains are produced and secreted. These cross-reactive antibodies may elicit the entry of viruses into macrophages in a Fc receptor-mediated manner, and consequently, viruses undergo rapid replication and release, resulting in immune dysregulation, and severe illness in patients with COVID-19 (52). The potential role of SARS-CoV-2-induced IgG antibodies in promoting neuroinflammation in SARS-CoV-2 infection should also be mentioned, as ADE occurrence may involve microglia cells following the binding of Abs to Fc receptors expressed on these cells.

As the mTOR pathway plays a fundamental role in B-cell development via the control of BCL6 expression in B cells from the germinal center (53), it is reasonable to believe that the use of inhibitors of mTOR, i.e., rapamycin and 'rapalogs', could reduce the populations of antigen-specific memory B cells and limit the occurrence of ADE in SARS-CoV-2 infected patients. This further strengthen the rational for using mTOR inhibitors in COVID-19, as previously discussed (20). Indeed, by using an anti-signature computational approach, our analysis showed that the mTOR inhibitor, sirolimus, may be a candidate drug to be used in COVID-19 patients, which is in line with data on the activation of the phosphoinositol 3-kinase (PI3K)/AKT/mTOR pathway in response to the infection with another coronavirus, MERS-CoV (54). Also, mTOR has been recognized as a key factor in regulating the replication of viruses (36,54-57), and in patients with H1N1 pneumonia, early treatment with corticosteroids in combination with rapamycin has been associated with improvement in multiple organ dysfunction, virus clearance, and shortened time in ventilators (58).

Hence, the use of mTOR inhibitors may have many-fold advantages on the course od SARS-CoV-2 infection, which could improve lung pathology, but also, the peripheral manifestations of the disease, including the CNS.

Interestingly, our data suggest potential reasons for the gender differences in COVID-19 susceptibility (2). Indeed, the prevalence in men is between 55 and 68% (59) and increased clinical severity and mortality has been reported (60). Certainly, female-specific hormonal factors can be involved. In this regard, it is notable that among the upregulated genes in BALF from COVID-19 patients, 22 (Log(q) value=-4) are regulated by AR (androgen receptor), while 15 are regulated by ESR1 (estrogen receptor 1) (Log(q) value=-1.8). AR is known to play a key role in both innate and adaptive immune responses (61,62), and ESR1 has been recognized as a regulator of interferon production and anti-viral responses (63). These observations may underly the different clinical response to SAR-CoV-2 infection in women and men. It is important to note that selective estrogen receptor modulators. such as toremifene, have already been proposed as potential drugs to treat coronavirus infections (64). These observations point to biological processes that may explain the lower female incidence and lethality of SARS-CoV-2 infection, offering candidate therapeutic options in patients suffering from COVID-19.

Finally, we have to acknowledge some of the limitations of the present study. First, the differentially expressed genes, that we have prioritized in our study, and the deconvolution analysis have been obtained from a really small cohort of patients, hence the data may be biased, due to the high degree of interindividual variability that characterize SARS-CoV-2 infection. Lung-specific gene expression profiles from homo-geneous COVID-19 patients will allow to better identify prognostic predictors and tailored therapeutic strategies. Second, the deconvolution analysis of the immune populations does not allow to assess the functionality of the immune cells and their actual involvement in COVID-19 pathology.