The noncoding and coding transcriptional landscape of the peripheral immune response in patients with COVID‐19

2.1. Study approval

The study was approved by the Ethics Committee of the Huo Shen Shan Hospital of Wuhan. All blood samples for multi‐transcriptome sequencing were existing samples that were collected during standard COVID‐19 treatment process, with no extra burden posed.

2.2. Sample collection

Whole blood was obtained from six severe and six moderate COVID‐19 patients at Huo Shen Shan Hospital of Wuhan during standard diagnostic tests. Sample collection criteria included age ≥18 years and admission to Huo Shen Shan Hospital (wards and ICU) with a positive result in SARS‐CoV‐2 nasopharyngeal swab RT‐PCR test. For controls, blood was collected from four uninfected adult donors with a negative nasopharyngeal swab. The classification of COVID‐19 severity was based on WHO guidelines. The classification of acute respiratory distress syndrome (ARDS) was based on the Berlin criteria (acute onset of hypoxemic respiratory failure with a PaO₂/FiO₂ < 300 on at least 5 cm of positive end‐expiratory pressure, bilateral infiltrates on chest X‐ray). Ages were shown as ranges to protect the privacy of patients and healthy controls. All donors were asked for consent for genetic research.

2.3. mRNA and miRNA isolation and purification

Erythrocytes were removed from human whole blood with Erythrocyte lysis buffer (Cat. No.: AR1118, BOSTER). Small RNAs (<200 nt) and large RNAs (>200 nt) were isolated using the NucleoSpin miRNA kit (Macherey‐Nagel, Düren, Germany) following the manufacturer's instructions. RNA purity was checked using the NanoPhotometer spectrophotometer (IMPLEN, CA). RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, CA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA).

2.4. Library construction and sequencing

For rRNA‐depleted RNA‐seq sample preparations, we used 3 μg RNA per sample. Sequencing libraries were generated using the rRNA‐depleted RNA by NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB) following manufacturer's instructions. The libraries were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired‐end reads were generated.

For miRNA‐seq, we used 1 μg RNA per sample for small RNA library construction. Sequencing libraries were generated using NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB) following manufacturer's instructions. To acquire sufficient sequencing coverage, four samples were combined into one lane and two technical replicates were ran for each library using multiplexing. T cluster generation was performed on a Flow Cell v3 (TruSeq SR Cluster Kit v3; Illumina) using cBOT. The library preparations were sequenced on an Illumina Hiseq 2500 platform and 50 bp single‐end reads were generated.

2.5. Data processing

Quality control processes consisted of adapter trimming, low‐quality reads removal with cutadapt software (version 2.10; https://cutadapt.readthedocs.io/en/stable/), and rRNA and tRNA removal. The total rRNA proportion indicated the quality of our samples with all samples (small and large RNA) containing rRNA less than 5%. All clean large RNA data were mapped to the human genome GRCh38 using HISAT2 (version 2.2.0; http://daehwankimlab.github.io/hisat2/). miRNA data were mapped to the miRBase (version 22; http://www.mirbase.org/) with Bowtie software (version 1.2.3; www.sourceforge.net/projects/bowtie-bio/files/bowtie), allowing 0 mismatch and mapping with the proximity of mature miRNAs. Bam files were sorted by Samtools (version 1.9; http://samtools.sourceforge.net/index.shtml). Gene counts were generated using the featureCounts program, part of the Subread package (version 2.0.0; http://subread.sourceforge.net/). miRNA counts were summarized with perl scripts written by authors of this paper.

The quality control of raw read data was done by FastQC (version 0.11.9) and multiQC (version v1.8). For rRNA‐depleted RNA‐seq data, the average number of reads per sample is about and the quality score for each sample is between 33 and 37, indicating good raw data quality. For miRNA‐seq data, two peaks at 22 nt and 33 nt were observed in read length distribution, with the first peak greater than the second peak. All RNA sequencing data exhibited rRNA alignment rates smaller than 5%. Sequencing data are available in the National Genomic Data Center (NGDC) (primary accession number HRA000238).

2.6. miRNA‐mRNA and lncRNA‐miRNA network construction

Both miRNA‐mRNA and lncRNA‐miRNA networks were constructed with Cytoscape (version 3.7.2; https://cytoscape.org/) based on their interaction. The R package “multimiR” (version 2.3; http://multimir.org) were used to identify miRNA‐mRNA interaction. We filtered every miRNA‐mRNA pair using the most stringent criteria, including validation by luciferase experiments and functional miRNA‐target interaction (MTI) tests. Database ENCORI (version; http://starbase.sysu.edu.cn/) were used to identify lncRNA‐miRNA interaction validated by cross‐linking and immunoprecipitation (CLIP)‐sEquation (≥5).

2.7. Single‐cell RNA‐seq computational pipelines and analysis

The R package Seurat (version 3.0; https://satijalab.org/seurat/) was used for single‐cell analysis including normalization, scaling, dimensionality reduction, clustering, transcriptome analysis, and visualization. Peripheral blood mononuclear cells (PBMCs) single‐cell RNA sequencing (scRNA‐seq) data (count matrix) from Wilk et al and nasopharynx scRNA‐seq data (Seurat object) from Chua et al were obtained from the provided available data sources for co‐analysis. PBMC data scaling, integration, clustering, and dimensionality reduction were performed following the R scripts provided by Wilk et al. Seurat's implementation of the Wilcoxon rank‐sum test (FindMarkers()) was used to determine differentially expressed genes (DEGs) for each cluster of both datasets. Cellular identity was determined based on marker gene expression presented in Supporting Information figures during re‐clustering.

2.8. Statistic

The power calculation was performed with the RNAseqPS tool ¹⁴ (https://cqs-vumc.shinyapps.io/rnaseqsamplesizeweb/). The number of patients offer sufficient power to detect twofold changes in gene expression based on depth and coverage of our sequencing data. The total quantity of each RNA was used to calculate the number of mapped reads per kilobase per million reads (RPKM). differentially expressed miRNAs (DEmiRs), DEGs, and differentially expressed lncRNAs (DElncRs) were reported using an adjusted P value threshold of .05, and a minimum fold change (FC) of 2. The adjusted P values were obtained from Independent Hypothesis Weighting as implemented in the DESeq2 package. The P‐values shown in each box plot figures are exact two‐sided by one‐way ANOVA with post hoc comparisons by Tukey's test using GraphPad Prism 8. Exact two‐sided P‐values and the 95%CI by Pearson correlation coefficients for each correlation are shown in scatter plots by GraphPad Prism 8.

3.1. Sample collection and transcriptome sequencing

RBC‐depleted whole blood samples were obtained from laboratory‐confirmed moderate (n = 6) and severe COVID‐19 patients (n = 6), as well as healthy controls (n = 4) (Figure 1A). Table 1 shows the detailed clinical characteristics of patients. Notably, the three groups exhibited no statistically significant difference in age between each other (Figure S2B).

We applied small RNA‐seq and rRNA‐depleted RNA‐seq on each blood sample. After quality control, reads were mapped to miRbase and GRCh38 genome, respectively. The count matrixes derived from the latter were further separated into protein‐coding mRNAs and lncRNAs.

We identified differentially expressed miRNAs (DEmiRs), DEGs, and DElncRs from comparisons between moderate‐healthy (M‐H), severe‐healthy (S‐H), and severe‐moderate (S‐M) (Figure S1). The list of DEmiRs, DEGs, and DElncRs was determined using adjusted P values (q‐value < 0.05) and FC ratios (|log2FC| ≥ 1) (Table S1‐S3).

Considering the total number of identified mRNAs and ncRNAs, miRNAs had the most evident alterations, supporting their high sensitivity as potential biomarkers (Figure S2A).

3.2. Identification of biomarker miRNAs by small RNA‐seq

The DEmiRs of each group are depicted in the heatmap (Figure S2B) and subdivided into different types based on their expression tendencies (Figure 1B). Four types of miRNAs were defined as candidate miRNA biomarkers (Figure 1C): (a) miRNAs consistently downregulated, including miR‐146a‐5p, miR‐21‐5p, and miR‐142‐3p (Figure 1D); (b) miRNAs consistently upregulated, including miR‐3605‐3p (Figure 1E); (c) miRNAs upregulated only in patients with severe COVID‐19 with no statistically significant difference in M‐H comparison, including miR‐15b‐5p, miR‐486‐3p, and miR‐486‐5p (Figure 1F); and (d) miRNAs downregulated only in severe cases, including miR‐181a‐2‐3p, miR‐31‐5p, and miR‐99a‐5p (Figure 1G).

Functional enrichment analysis of the predicted target genes for several representative biomarker miRNAs showed high correlation with inflammation and antiviral immune responses (Figure 1H). Processes including “virus binding,” “virus process,” and “defense response to virus” implied miRNA engagement in viral infection. These miRNAs were also associated with several Toll‐like receptor (TLR) signaling pathways, as well as production of and response to interferon.

3.3. Bulk RNA‐seq suggests biomarker miRNAs as potential therapeutic targets

To better understand the function of biomarker miRNAs, we retrieved their validated downstream mRNAs, analyzed the miRNA‐mRNA correlation, and constructed an integrated miRNA‐mRNA regulatory network with “multimiR” package (Figure 2A). Of note, we filtered every miRNA‐mRNA pair with the most stringent criteria, including validation by luciferase experiments and functional MTI tests.

Consistent downregulation of miR‐146a‐5p, miR‐21‐5p, and miR‐142‐3p promotes inflammatoty process. ¹⁵ , ¹⁶ , ¹⁷ The miR‐146a‐5p negatively correlated with downstream target mRNAs interleukin‐1 receptor‐associated kinase 1 (IRAK1), IRAK2, and tumor necrosis factor receptor‐ associated factor 6 (TRAF6), which participate in the nuclear factor‐κB (NF‐κB) pro‐inflammatory signaling pathway ¹⁸ , ¹⁹ , ²⁰ (Figure 2B). miR‐21‐5p may directly target IRAK1 and chemokine C‐C motif chemokine ligand 20 (CCL20), which was upregulated in inflamed airway epithelium ²¹ (Figure 2C). Decreased miR‐142‐3p induces production of glycoprotein 130 (gp130), an activator of JAK/STAT signaling pathway, by binding to interleukin 6 signal transducer (IL6ST) mRNA ²² (Figure 2D).

miRNAs that were up‐ or downregulated only in severe cases may contribute to COVID‐19 deterioration. Upregulated miR‐15b‐5p seemed to play dual roles. First, it negatively correlated with IFNG and CD69 that were involved in T‐cell function and activation ²³ , ²⁴ (Figure 2E). Moreover, miRNAs can promote RNA virus replication by binding to and stabilizing the viral genome. A recent research identified miR‐15b‐5p as the most likely candidate to target SARS‐CoV‐2 genome with the highest target score and binding sites ²⁵ (Figure S2C). Upregulation of miR‐15b‐5p could accelerate intracellular viral replication, promote cell‐to‐cell dissemination, mediate virus‐induced transcriptome changes, and ultimately intensify the severity of COVID‐19.

Other upregulated candidate miRNAs also contributed to COVID‐19 pathogenesis. miR‐486‐5p not only targets neuropilin 2 (NRP2) encoding inflammatory inhibitor neuropilin‐2, but also represses OTUD7B to induce excessive inflammation in lung ²⁶ , ²⁷ , ²⁸ (Figure 2F). Another miRNA miR‐486‐3p directly targets MAF, and downregulation of MAF may result in immune response dysregulation ²⁹ , ³⁰ (Figure 2G).

Among the miRNAs downregulated in severe cases was miR‐181a‐2‐3p, a serum biomarker of chronic obstructive pulmonary disease. ³¹ Downregulation of miR‐181a‐2‐3p is also associated with enhanced TLR4 and CXCL8 expression ³² , ³³ (Figure 2H). Another downregulated miRNA was miR‐99a‐5p targeting the proinflammatory genes, insulin like growth factor 1 receptor (IGF1R) ³⁴ , ³⁵ and myotubularin‐related protein 3 (MTMR3), which were reported to induce weaker antiviral immunity ³⁶ , ³⁷ (Figure 2I).

Overall, the results of our detailed functional miRNA‐mRNA analysis suggest that the candidate biomarker miRNAs identified may contribute to COVID‐19 pathogenesis and serve as therapeutic targets.

3.4. Validation of miRNAs as therapeutic targets by single cell RNA‐seq

To further validate the function of candidate miRNAs, we used the scRNA‐seq data from Wilk et al, pertaining to the sequencing of PBMCs from severe cases and healthy controls. ³ The samples from Wilk et al study could result in different DEGs with ours because of the absence of neutrophils. We successfully identified 13 major peripheral blood cell clusters (Figure 3A). Monocytes appeared to be the most remodeled cells (Figure S3A). The major groups of immune cells were re‐clustered for further analysis (Figure S3B).

We first focused on consistently downregulated miRNAs. Our analysis revealed broad upregulation of STAT1 targeted by miR‐146a‐5p that encodes a key element of the JAK/STAT pathway ³⁸ , ³⁹ (Figure 3B). As for miR‐21‐5p, MYC, a T‐cell activation marker, displayed upregulation specifically in naïve T cells ⁴⁰ , ⁴¹ (Figure 3C and D). IL6ST, a target of miR‐142‐3p, was upregulated mainly in naïve T cells, suggesting increased sensitivity to interleukin (IL)‐6 signaling ⁴² (Figure 3E). Upregulation of STAT1, MYC, and IL6ST was also observed in our bulk RNA‐seq data with smaller fold changes and larger adjusted P value (Figure S3C).

It was improper to use the scRNA‐seq data by Wilk et al. that only addressed the S‐H comparison to validate the function of miRNAs differentially expressed in the S‐M comparison. Neutrophils were also excluded from PBMC samples. Thus, we used the scRNA‐seq data of Chua et al obtained in cells of the nasopharyngeal area of severe and moderate COVID‐19 patients, and healthy controls. ⁴ We focused on cells transported from the peripheral blood to the nasopharyngeal area, whose transcriptome tended to be regulated by blood biomarker miRNAs (Figure 3F). The targets of miR‐15b‐5p, IFNG, and CD69 were downregulated in nasopharyngeal CD8+ T cells in S‐M comparison. We also observed the upregulation of miR‐99a‐5p target, MTMR3 as well as miR‐181a‐2‐3p targets, TLR4 and CXCL8, in neutrophils in the context of severe disease (Figure S3D and E).

Re‐examining miR‐146a‐5p, miR‐21‐5p, and miR‐142‐3p targets with Chua et al data successfully identified upregulation of IRAK1 and IRAK2 along with other targets mainly expressed in neutrophils ⁴³ (Figure S4).

Overall, the combined analysis of our miRNA data and two published scRNA‐seq datasets further validated these candidate biomarker miRNAs as potential therapeutic targets in specific cell types.

3.5. lncRNAs could be the upstream “sponges” inhibiting miRNA function

lncRNAs may act as miRNA sponges that bind to specific miRNA sites, reduce miRNA‐mRNA interaction, and inhibit the regulatory function of miRNAs. ⁸ To provide a comprehensive overview of the upstream regulator of our potential therapeutic target miRNAs, we established a lncRNA‐miRNA network based on lncRNA‐miRNA CLIP results.

DElncRs in each group were shown in the heatmap (Figure S5A). We explored lncRNA‐miRNA interactions using the StarBase database, set up filter criteria (CLIP ≥ 5), and applied it to our lncRNA data (Figure 4A). Nuclear paraspeckle assembly transcript 1 (NEAT1) seemed to dominate the networks with four associated miRNAs, implying its critical role in COVID‐19 pathogenesis. Correspondingly, NEAT1 showed increased expression across S‐M and S‐H comparisons, with no significant change in the M‐H comparison (Figure 4B). Moreover, NEAT1 might increase the production of inflammatory cytokines, IL‐6 and CXCL8, which confirmed its critical role in inflammation ⁴⁴ , ⁴⁵ (Figure 4C).

We identified several lncRNA‐miRNA pairs with opposite expression tendencies and biological functions. Both miR‐146a‐5p and miR‐142‐3p were negatively correlated with a canonical inflammatory inhibitor, metastasis associated lung adenocarcinoma transcript 1 (MALAT1) ⁴⁶ (Figure 4D). MALAT1 could absorb miR‐146a‐5p and miR‐142‐3p to repress their anti‐inflammatory function. ⁴⁷ , ⁴⁸ miR‐15b‐5p could be targeted by multiple lncRNAs, among which inflammatory inhibitors, long intergenic nonprotein coding RNA 649 (LINC00649) and small nucleolar RNA host gene 1 (SNHG1), showed opposite expression patterns ⁴⁹ (Figure 4E).

Targeting these lncRNAs could affect the silencing function of their downstream miRNAs, which provides another treatment strategy for COVID‐19.

3.6. Cytokines, TLR, and marker gene profiles

Elevated inflammatory mediators play a crucial role in fatal pneumonia caused by SARS‐CoV‐2. ² To better understand the potential effect of each type of cytokines, we sorted genes into six categories including “interleukin,” “interferon,” “chemokine,” “tumor necrosis factor,” “other,” and “receptor” of which “interleukin” and “chemokine” displayed the most evident changes in gene expression (Figure 5A).

ILs exhibited the largest fold changes (Figure 5B). Consistent with previous large‐scale clinical studies, we observed significant upregulation of IL10 and IL6 in COVID‐19 patients, suggesting the ongoing hyperactivation of inflammatory effects. ⁵⁰ , ⁵¹ , ⁵² In addition, we observed a consistent downregulation of IL4 that limits tissue damage during the immune response. ⁵³

We detected consistent upregulation of CXCL2 critical for recruiting neutrophils into inflamed lungs. ⁵⁴ The expression of another important proinflammatory cytokine CCL20 was upregulated. We also observed increased expression of CXCL1 and CXCL6 in S‐M comparison, indicating an excessive activation of inflammatory effects in severe cases ⁵⁵ , ⁵⁶ (Figure S6).

The downregulation of IFNG and CD69 in severe cases (probably regulated by miR‐15b‐5p) led us to T‐cell exhaustion analyses. This downregulation was also observed in PBMCs, bronchoalveolar lavage fluid (BALF), and nasopharynx cells from severe COVID‐19 patients. ² , ⁴ , ⁵⁷ T‐cell exhaustion‐associated genes were separated into “Activation,” “Function,” and “Regulation” categories (Figure 5C). CD69, IFNG, TNF, IL2, TLR4, GZMB, and PRF1 were downregulated in severe COVID‐19 patients compared to moderate cases, suggesting loss of T‐cell activation, cytokine secretion, and cytotoxicity. ⁵⁸ , ⁵⁹ Moreover, increased expression of IL10 may induce T‐cell dysfunction. ⁵⁸

Analyses of immune cell markers reflected decrease of T‐cell numbers in severe cases with downregulated CD3, CD4, and CD8, which was corroborated by large‐scale clinical studies ⁵⁰ (Figure 5D). We also observed a significant upregulation of MZB1, which agrees with a recent study demonstrating the unique role of plasmablasts in severe COVID‐19.

Because our functional enrichment analysis revealed an enrichment in miRNAs involved in TLR signaling, we set out to investigate the expression of TLR genes (Figure 5E). Most TLR mRNAs, including TLR4, TLR5, and TLR8 were upregulated in the severe group, consistent with viral and bacterial co‐infection that was common among severe cases. Interestingly, two antiviral immunity‐associated TLRs, TLR3 and TLR7, showed opposite expression profiles, which may account for diminished viral clearance and disease aggravation in severe cases. ⁶⁰

3.7. Global analysis of DEGs during COVID‐19 by bulk RNA‐seq

Compared to miRNAs, the expression of mRNAs was more complex and characterized by less distinct DEG clusters (Figure S5B). Similar to our miRNA analysis, we searched for DEGs across the three comparisons (M‐H, S‐H, and S‐M) and obtained a list of DEGs possibly associated with disease pathogenesis (Figure 6A). Four groups of DEGs were identified (Figure 6B): (a) genes consistently upregulated in COVID‐19 patients, including the marker of cell proliferation, marker of proliferation Ki‐67 (MKI67), ⁶¹ and regulator of macrophage function, forkhead box M 1 (FOXM1) ⁶² (Figure 6C); (b) genes consistently downregulated, including methyltransferase like 21C (METTL21C) regulating NF‐κB signaling and the expression of IL10 ⁶³ (Figure 6D); (c) genes exclusively upregulated in severe COVID‐19 patients, including absent in melanoma 2 (AIM2), PIM1, and angiotensin I converting enzyme 2 (ACE2) (Figure 6E); cytoplasmic DNA can elicit AIM2 expression to induce cell death ⁶⁴ while PIM1 kinase promotes airway inflammation; ⁶⁵ ACE2 will be discussed in detail later; (d) genes exclusively downregulated in severe COVID‐19 patients, including CD8+ T‐cell inhibitor CD248 ⁶⁶ and cytotoxicity marker, granulysin (GNLY) ⁶⁷ (Figure 6F).

Functional enrichment analysis identified significant enrichments in the “response to bacterium” reflecting bacterial co‐infections and processes related to lymphocyte activation, proliferation, and regulation (Figure 6G). KEGG pathway analysis identified the “cytokine and cytokine receptor interaction” and “viral protein interaction with cytokine and cytokine receptor,” suggesting hypercytokinemia caused by severe infection (Figure 6H). Upregulation of “NF−κB signaling pathway” indicated the activation of common upstream pathways regulating cytokine production. The upregulated “IL‐17 signaling pathway” agreed with elevated levels of Th17 cells observed in COVID‐19 patients. ⁶⁸

Notably, ACE2 was only upregulated in S‐M and S‐H comparisons (Figure 6E). ACE2 encodes the angiotensin I converting enzyme 2, a cell receptor considered vital for the entry of SARS‐CoV‐2 in host cells. ⁶⁹ The scientific community has proposed that ACE2 expression can be induced by interferon via STAT1 signaling ⁴ and correspondingly, we observed the upregulation of IFNL2 (Figure 6I). Interestingly, the fact that ACE2 was mainly upregulated in epithelial cells coincided with the exclusive expression of interferon (IFN)‐λ receptor at the surface of epithelial cells. ⁷⁰ This observation may imply that upregulation of ACE2 is induced by IFN‐λ and mediated by STAT1 signaling in lung tissues.