Plasma cell-free RNA characteristics in COVID-19 patients.

The pathogenesis of COVID-19 is still elusive, which impedes disease progression prediction, differential diagnosis, and targeted therapy. Plasma cell-free RNAs (cfRNAs) carry unique information from human tissue and thus could point to resourceful solutions for pathogenesis and host-pathogen interactions. Here, we performed a comparative analysis of cfRNA profiles between COVID-19 patients and healthy donors using serial plasma. Analyses of the cfRNA landscape, potential gene regulatory mechanisms, dynamic changes in tRNA pools upon infection, and microbial communities were performed. A total of 380 cfRNA molecules were up-regulated in all COVID-19 patients, of which seven could serve as potential biomarkers (AUC > 0.85) with great sensitivity and specificity. Antiviral (NFKB1A, IFITM3, and IFI27) and neutrophil activation (S100A8, CD68, and CD63)-related genes exhibited decreased expression levels during treatment in COVID-19 patients, which is in accordance with the dynamically enhanced inflammatory response in COVID-19 patients. Noncoding RNAs, including some microRNAs (let 7 family) and long noncoding RNAs (GJA9-MYCBP) targeting interleukin (IL6/IL6R), were differentially expressed between COVID-19 patients and healthy donors, which accounts for the potential core mechanism of cytokine storm syndromes; the tRNA pools change significantly between the COVID-19 and healthy group, leading to the accumulation of SARS-CoV-2 biased codons, which facilitate SARS-CoV-2 replication. Finally, several pneumonia-related microorganisms were detected in the plasma of COVID-19 patients, raising the possibility of simultaneously monitoring immune response regulation and microbial communities using cfRNA analysis. This study fills the knowledge gap in the plasma cfRNA landscape of COVID-19 patients and offers insight into the potential mechanisms of cfRNAs to explain COVID-19 pathogenesis.

COVID-19 is currently threatening global health. As of November 2021, there are more than 253 million COVID-19 patients and 5.1 million deaths reported, and case numbers are still rising in many countries (WHO Coronavirus [COVID-19] Dashboard; https://covid19.who.int/). There remains an urgent need for elucidating the pathogenesis of COVID-19. Recent studies suggested that the cytokine storm syndromes contribute to the high mortality rate of COVID-19, and interleukin 6 (IL6)/interleukin 6 receptor (IL6R) signaling plays a crucial role in the cytokine storm (Scheller and Rose-John 2006; Guan et al. 2020). SARS-CoV-2 infection activates the function of immune cells and releases IL6, therefore causing an amplification cascade of cytokine release syndromes (CRSs) and finally resulting in the pathophysiology of COVID-19 (Moore and June 2020). Down-regulated miR-451a as a feature of the plasma cfRNA landscape reveals regulatory networks of IL6/IL6R-associated cytokine storms in COVID-19 patients (Yang et al. 2021).

Cell-free circulating RNAs (cfRNAs) consist of a variety of RNA molecules, such as messenger RNAs (mRNAs) (Ng et al. 2002), microRNAs (miRNAs) (Pritchard et al. 2012) and long noncoding RNAs (lncRNAs) (Kumarswamy et al. 2014). cfRNAs present in various body fluids, such as blood, serum, and plasma, and are involved in multiple physiological and pathological processes (Pös et al. 2018). Although the mechanisms of how cfRNAs are released to the bloodstream remain unclear, they can serve as valuable biomarkers to indicate dynamic changes of different tissues and enable the detection of specific genetic alterations (Tzimagiorgis et al. 2011). Compared with cell and tissue samples, plasma cfRNAs carry unique information from the human tissue or pathologic area, including pathological states, disease biomarkers (Larson et al. 2021), gene regulatory networks, and microbial co-infections, which can provide resourceful solutions for pathogenesis. Due to their typical advantages of real-time and comprehensive monitoring, cfRNAs have the potential to change the ways of monitoring and predicting patients’ status.

Transfer RNAs (tRNAs) serve as the carrier of amino acids and participate in protein synthesis. tRNAs can recognize codons on mRNA by its anticodon. Within the coding sequence of mRNAs of a certain organism, it is observed that some synonymous codons are used more frequently than other synonymous codons, and this phenomenon is called codon bias (Grantham et al. 1980). Genes with high expression levels tend to use optimal codons to improve translation accuracy and efficiency, showing high codon usage preference (Hanson and Coller 2018). Moreover, codon bias can also affect protein folding and differential regulation of protein expression (Mitra et al. 2016). For example, in humans, codons with G or C at the third base position can stabilize mRNA (Hia et al. 2019). In contrast, SARS-CoV-2 preferred fewer G or C nucleotides at the third position and this codon preference may disturb protein synthesis of the host cells (Hou 2020). The tRNA pool also changes dynamically with a spatial-temporal expression pattern (Torrent et al. 2018; Torres et al. 2019; Wei et al. 2019). For example, during influenza A virus (FLUAV) infection, the expression levels of codons preferred by FLUAV were increased, which may promote the replication and translation of viral genes (Pavon-Eternod et al. 2013).

Microbial co-infections (viruses, bacteria, and fungi) often occur in respiratory infections (Langford et al. 2020), and this also happens in COVID-19 (Hughes et al. 2020; Massey et al. 2020; Shen et al. 2020; Zhu et al. 2020b). Previous studies found that common co-infective pathogens in COVID-19 patients include Klebsiella pneumoniae (Massey et al. 2020; Zhu et al. 2020b), Haemophilus influenzae (Zhu et al. 2020b), Pseudomonas aeruginosa (Hughes et al. 2020), and Moraxella catarrhalis (Massey et al. 2020) which raise the difficulty of diagnosis, treatment, and prognosis of COVID-19 and even aggravate the symptoms (Shen et al. 2020; Zhu et al. 2020b). The analysis of cfRNA is a noninvasive approach for molecular diagnostics (Pös et al. 2018); currently no evidence of microbial co-infections derived from cfRNA analysis has been reported that impeded the integrated evaluation of disease status.

At present, advanced omics research has revealed different perspectives of COVID-19 (Alonso and Diambra 2020; Hughes et al. 2020; Massey et al. 2020; Shen et al. 2020; Zhu et al. 2020a). For instance, one study unveiled characteristic protein and metabolite changes in the sera of severe COVID-19 patients (Shen et al. 2020), and another revealed the unique immune response pathways in COVID-19 patients compared to influenza patients through single-cell RNA sequencing of peripheral blood mononuclear cells (Zhu et al. 2020a). However, the identification of COVID-19 biomarkers and analysis of the disease process from the perspective of cfRNA remain to be investigated. Here, we performed PALM-seq to provide a plasma cfRNA landscape and open up new perspectives for understanding COVID-19 pathogenesis.

Results

Characteristics of cfRNA landscape in COVID-19 patients

Serial plasma was collected from eight healthy donors and 37 COVID-19 patients (19 mild and 18 severe patients) and used for cfRNAs sequencing by PALM-seq (Fig. 1A; Supplemental Table S1; Yang et al. 2019). Based on the clinical symptoms, two distinct subgroups of COVID-19 patients were enrolled, including the mild and severe groups. Patients with severe pneumonia who were admitted to the ICU and required mechanical ventilation were enrolled in the severe group; patients with a mild clinical presentation (mainly with fever, cough, malaise, and headache) were enrolled in the mild group. All available COVID-19 patients were enrolled in this study, and healthy donors aged between 20∼60 were selected as controls. We first aligned the reads to the SARS-CoV-2 genome to determine the presence of SARS-CoV-2 RNAs in plasma samples. On average, <0.008% of the reads could be mapped to SARS-CoV-2 (Supplemental Fig. S1A), and most of the reads mapped to SARS-CoV-2 (31.8%) were derived from four severe patients, which is consistent with recent studies reporting that SARS-CoV-2 RNA was rarely detected in plasma except certain severe patients (Ling et al. 2020; Wölfel et al. 2020; Yang et al. 2021).

Figure 1.

The plasma cfRNA landscape reveals pathogenesis of COVID-19 patients. (A) Schematic diagram showing samples and analysis procedure of this study. The cfRNAs were collected from eight healthy donors, 19 mild COVID-19 patients, and 18 severe COVID-19 patients. (B) Bar plot and pie plot showing the composition of detected genes in all samples. Data are shown as mean + SD (n = 156). (C) Heat map showing fold changes of 380 up-regulated cfRNA genes and 117 down-regulated cfRNA genes in COVID-19 patients. Fold changes of the 380 up-regulated genes were relative to healthy donors, and those of the 117 down-regulated genes were relative to COVID-19 patients. Blue and red represent log2-transformed fold changes < 0 and > 0, respectively. (D) Overlapping of up-regulated genes in mild and severe COVID-19 patients. (E) GO enrichment analyses of 231 up-regulated cfRNA genes in COVID-19 patients and 137 cfRNA genes specific up-regulated in severe COVID-19 patients. (F) Box plot showing expression levels (TPMs) of NFKBIA and S100A8 in healthy donors, mild, and severe COVID-19 patients. (G) Box plot showing the neutrophil count in mild and severe COVID-19 patients. Asterisks indicate statistically significant differences; (***) P < 0.001.

Next, all the reads were aligned to the human genome and 27,738 expressed genes were analyzed, including 71.2% protein-coding genes, 6.1% miRNA genes, 1.3% lncRNA genes, and 21.4% other noncoding genes (Fig. 1B). By comparing the data of COVID-19 patients and healthy donors, 380 up-regulated genes were identified, including 231 genes up-regulated in both mild and severe patients, 12 and 137 genes up-regulated in mild or severe patients, respectively (Fig. 1C,D; Supplemental Table S2). The hierarchical clustering results showed that these 380 up-regulated genes and 117 down-regulated genes could distinguish COVID-19 patients from healthy donors (Fig. 1C; Supplemental Table S2). We also calculated the receiver operating characteristic (ROC) curve of each up-regulated gene and found that seven of them exhibited area under the curve (AUC) > 0.85 (Bowers and Zhou 2019), indicating that these genes could be served as potential biomarkers for prediction of disease. (Supplemental Table S3). In addition, we performed logistic regression on the 380 up-regulated genes and screened 57 genes (with 98.1% accuracy) as the optimal gene set for COVID-19 prediction (Supplemental Table S3). Consistent with ROC analysis, we found that all of these seven biomarkers were included in this gene set. These results further illustrate that these biomarkers were significantly differentially expressed between COVID-19 patients and healthy donors. Gene Ontology (GO) function annotation revealed that the up-regulated 231 genes in both mild and severe patients were mainly enriched in anti-virus-related pathways (Fig. 1D,E), indicating that anti-virus response was activated in both mild and severe patients. Representative genes such as NFKB1A, IFITM3, and IFI27, were all important regulators of anti-virus processes (Fig. 1F; Supplemental Fig. S1B; Hayden et al. 2006; Tang et al. 2017; Kenney et al. 2019; Spence et al. 2019). Previous studies have also found that IFI27 was up-regulated in patients with COVID-19 (Yang et al. 2021), which is consistent with our results. The functions of the 137 genes that were specifically up-regulated in severe patients were mainly enriched in neutrophil activation–related processes (Fig. 1D,E), potentially indicating a strong inflammatory response in severe patients, which is consistent with recent studies reporting that disease severity–specific neutrophil signatures in blood stratify COVID-19 patients (Weiland et al. 1986; Aschenbrenner et al. 2021). These genes include well-known key genes for neutrophil activation such as S100A8, CD68, and CD63 (Fig. 1F; Supplemental Fig. S1B; Skubitz et al. 2000; Ryckman et al. 2003; Amanzada et al. 2013; Wang et al. 2018). A previous study found that S100A8, CD68, and CD63 had low transcription start site coverage in severe COVID-19 patients compared with mild patients, indicating that these genes had higher expression levels in severe COVID-19 patients (Chen et al. 2021). Measurement of neutrophil count was performed in COVID-19 patients, and as expected, it was significantly higher in severe patients than in mild patients (Fig. 1G).

To study the transcriptional dynamics of cfRNAs during the disease process, the whole clinical course from admission to discharge was divided into three stages (early, middle, and late) (also see Supplemental Table S1). Time-course analysis was performed, and the results showed that 682 and 431 genes exhibited four and three expression patterns in severe and mild patients, respectively (Supplemental Fig. S1C,D). Cluster 1 and 2 genes of severe (S-1 and S-2) and mild patients (M-1 and M-2) exhibited similar patterns whose expression levels were increased at the middle stage and then decreased at the late stage. GO enrichment analyses revealed that these genes were mainly related to neutrophil degranulation (Supplemental Fig. S1E), further supporting the clinical finding that inflammation responses were reduced at the late stage (Jaovisidha et al. 1999). Cluster 4 genes in severe patients (S-4) were gradually down-regulated from early to late stage, and their functions were enriched in mRNA splice-related pathways (Supplemental Fig. S1E), suggesting that the virus infection may cause changes in RNA splice, affecting the processing and maturation of host mRNA (Dubois et al. 2014; Chauhan et al. 2019). Cluster 3 genes in mild patients (M-3) were enriched in the type I interferon signaling pathway (Supplemental Fig. S1E), indicating that the inflammatory response was gradually reduced.

Noncoding regulatory RNAs may aggravate IL6- and IL6R-induced cytokine storm in COVID-19 patients

The IL6/IL6R signaling cascade plays an important role in the elevated inflammatory responses (cytokine storm) of COVID-19 patients (Scheller and Rose-John 2006; Guan et al. 2020). Given that inflammatory responses were hyperactivated in COVID-19 patients (Fig. 1; Supplemental Fig. S1), we next investigated the potential regulatory mechanism. First, we examined the mRNA levels of IL6/IL6R and found no significant differences between COVID-19 patients and healthy donors (Fig. 2A). In contrast, increased plasma IL6 protein concentrations were observed in most COVID-19 patients (78%, 14 of 18) (Fig. 2A), in agreement with a previous study (Shen et al. 2020). These results raised the possibility that IL6/IL6R function might be enhanced at the translational level.

Figure 2.

Noncoding regulatory RNA responses of COVID-19 patients to SARS-CoV-2 infection. (A) The mRNA expression levels (TPMs) of IL6 (left) and IL6R (middle), and the protein levels of IL6 (right) in COVID-19 patients’ plasma. The red line indicates the normal range of IL6 (0–7 pg/mL). (B) Predicted (left, dashed lines) and experimentally validated (right, solid lines) miRNA-target interactions. (C) Box plot showing the expression levels of hsa-let-7a-5p and hsa-let-7g-5p. (D) Bar plot showing the summation of the entire let-7 family expression levels. (E) Base-pairing interaction between let-7 family and IL6/IL6R (top) or GJA9-MYCBP (bottom). Target sites and seed sequences are highlighted in red. (F) Box plot showing the expression levels of GJA9-MYCBP. (G) Putative regulatory network of GJA9-MYCBP and let-7 family. Asterisks indicate statistically significant differences; (***) P < 0.001, (n.s.) not significant.

Noncoding RNAs (ncRNAs) play important roles in apoptosis, proliferation, cancer, and immune response (Zhang et al. 2012; Kopp and Mendell 2018). In animals, miRNAs can bind to the 3′ UTR of mRNAs and inhibit their translation (Wang et al. 2013). To investigate their roles during SARS-CoV-2 infection, we examined their expression levels and found that, compared with healthy donors, the global miRNA expression was down-regulated in COVID-19 patients, including 41 miRNAs that were down-regulated in both mild and severe COVID-19 patients (Supplemental Fig. S2A,B). To understand their potential functions, we predicted their target genes and found that several of them directly target IL6 and/or IL6R (Fig. 2B), both of which encode key drivers of the cytokine storm induced in COVID-19 (Scheller and Rose-John 2006; Iliopoulos et al. 2009; Guan et al. 2020; Moore and June 2020). Previous studies have experimentally validated that hsa-let-7a-5p and hsa-let-7g-5p could target IL6 mRNA and inhibit its translation (Iliopoulos et al. 2009; Huang et al. 2017). Consistent with this, we found that the expression levels of hsa-let-7a-5p and hsa-let-7g-5p, and also the entire let-7 family, were down-regulated in both mild and severe COVID-19 patients (Fig. 2C,D; Supplemental Fig. S2C). Moreover, hsa-let-7a-5p can also target IL6R (Sui et al. 2019), indicating that the let-7 family can inhibit the translation of both IL6 and IL6R (Fig. 2E). In addition to let-7 family members, we also found that hsa-miR-451a, hsa-miR-23a-3p, and hsa-miR-23b-3p were down-regulated in COVID-19 patients (Supplemental Fig. S2D), which has been shown to target IL6R (Fig. 2B; Supplemental Fig. S2E) and reduce its expression level (Aghaee-Bakhtiari et al. 2015; Liu et al. 2016). Previous studies have also found that hsa-miR-451a was down-regulated in patients with COVID-19 (Yang et al. 2021), which reconfirmed our results. Taken together, these results suggested that the down-regulation of these miRNAs may lead to excessive accumulation of the protein products of their target genes, thereby enhancing the cytokine storm caused by the IL6/IL6R cascade under SARS-CoV-2 infection (Fig. 2A,B).

Long noncoding RNAs can function as miRNA sponges by competing with their targets to bind miRNAs, which negatively regulate miRNA functions (Wang et al. 2013). In our data, we detected 4324 annotated lncRNA genes, including 1665 (38.5%) transcribed from the intergenic region and 2659 (61.5%) transcribed from the genic region (from both the same direction and the reverse direction overlapped with the annotated genes) (Supplemental Table S4). Among them, 1486 (34.4%) were fully covered by reads. We found lncRNA GJA9-MYCBP was up-regulated in severe COVID-19 patients (Fig. 2F), which carries a binding site of let-7 family members (Fig. 2E,G). Therefore, GJA9-MYCBP could bind to let-7 family members in COVID-19 patients (Fig. 2G) and increase the protein levels of IL6 and IL6R, and thereby promote the cytokine storm caused by IL6/IL6R.

tRNA pool perturbation toward SARS-CoV-2 biased codons in COVID-19 patients

In addition to translational efficiency, protein synthesis is also determined by tRNA pools which are dynamically regulated under different physiological conditions (Torrent et al. 2018; Torres et al. 2019; Wei et al. 2019), especially in RNA virus infections (Pavon-Eternod et al. 2013; Nunes et al. 2020). tRNA abundance affects the translational speed of both host and viral proteins, as the translation of viral proteins requires the host tRNAs (Pavon-Eternod et al. 2013; Lyons et al. 2018; Chen et al. 2020a). Here, 380 up-regulated genes have been identified in COVID-19 patients (Fig. 1C,D). We then investigated their compositions and found that, unexpectedly, ∼36% (138 of 380) were tRNA genes (Fig. 3A), representing 34% (138 of 404) of expressed tRNA genes (Fig. 3B). In contrast, only 1% of protein-coding genes (181 of 15,319) were up-regulated (Fig. 3B). This observation indicated that the global tRNA profiles were significantly altered in COVID-19 patients compared with healthy donors. We then compared the tRNA profiles between mild and severe COVID-19 patients and found that, among the 138 up-regulated tRNA genes, ∼95% (131 of 138) were up-regulated in both mild and severe COVID-19 patients (Supplemental Fig. S3A). The ROC curve further showed an AUC of 0.87 (Supplemental Fig. S3B), suggesting these 138 up-regulated tRNAs could differentiate COVID-19 patients from healthy donors. In addition, we also observed a high correlation (r = 0.90, P-value → 0) of all expressed tRNA genes between mild and severe COVID-19 patients (Fig. 3C).

Figure 3.

Aberrant tRNA pools and codon abundance in COVID-19 patients. (A) Pie plot showing the composition of the 380 up-regulated genes in COVID-19 patients compared with healthy donors. (B) Bar plot showing the proportion of up-regulated genes and codons in COVID-19 patients compared with healthy donors. (C) Scatterplot showing up-regulated (red) and down-regulated (blue) tRNA genes of COVID-19 patients relative to healthy donors. (D) Scatterplot showing up-regulated and down-regulated codons of COVID-19 patients relative to healthy donors. (E) PCA plot of codon abundance in all samples. The abundance of 47 detected codons shown in D was used to draw this plot. (F) Relative codon usage frequency of 19 up-regulated codons in up-regulated genes (left) and the SARS-CoV-2 genome (right). The codon usage frequency of the 19 up-regulated codons in other expressed genes (left) and human genome (right) was set to 1. (G) Nucleotide composition of 14 up-regulated codons which were more frequently used in the SARS-CoV-2 genome compared to the human genome. (H) Nucleotide composition of 11 up-regulated codons which were more frequently used in the up-regulated genes compared to the other expressed genes.

tRNAs with distinct anticodons participate in protein synthesis by carrying amino acids (Lyons et al. 2018); therefore, changes in tRNA expressions can lead to changes in anticodon abundance. To better understand its potential effects on protein synthesis, we merged the expression levels of tRNAs with the same anticodon and converted the anticodon sequences to the corresponding codon sequences. A total of 47 codons were identified and their global expression profiles exhibited a higher correlation (r = 0.95, P-value → 0) between mild and severe COVID-19 patients (Fig. 3D). Among them, 26 codons were up-regulated in COVID-19 patients, accounting for 55% of detected codons (Fig. 3B; Supplemental Fig. S3C). This result further supported that, compared with healthy donors, both the tRNA abundance and the corresponding codon abundance were significantly changed in COVID-19 patients. As a result, principal component analysis (PCA) using expression levels of all detected codons clearly separated COVID-19 patients from healthy donors with no obvious differences between mild and severe patients (Fig. 3E), which was in agreement with the correlation analysis (Fig. 3D). Therefore, we focused on the 26 up-regulated codons and further investigated their functions.

Codon abundance is affected by the codon usage frequency of highly expressed proteins (Goñi et al. 2012; Diambra 2017). Under SARS-CoV-2 infection, there are two types of highly expressed proteins: one is viral proteins, and the other is up-regulated host proteins, both of which may lead to the up-regulation of the 26 codons in COVID-19 patients. Thus, we calculated the fold change of codon usage frequency in the SARS-CoV-2 genome relative to the human genome and also, that in up-regulated genes relative to other expressed genes in COVID-19 patients, respectively (Fig. 3F). Of the 26 up-regulated codons, 19 were related to differences in codon usage frequency, including 14 codons (6 + 8) that were more frequently used in the SARS-CoV-2 genome compared to the human genome, and 11 codons (6 + 5) that were more frequently used in the up-regulated genes compared to the other expressed genes in COVID-19 patients (Fig. 3F). Although six codons were more frequently used in both conditions, the fold changes of SARS-CoV-2 genome were significantly higher than those of up-regulated genes (Fig. 3F). It is reported that, when the virus replicates in the host cell, a large number of the host tRNAs are up-regulated to increase the abundance of virus-preferred codons (Chen et al. 2020a). A similar phenomenon also happened here: the up-regulation of the 14 codons with higher codon usage frequency in the SARS-CoV-2 genome may facilitate its replication.

The increased expression of SARS-CoV-2-preferred AT3 codons may reduce the stability of the host mRNAs

Previous studies have found that codon usage bias also affects mRNA stability (Hia et al. 2019). Codons can be divided into two groups based on the third base: AT3 (A or T) or GC3 (G or C), and the bias of the third base is characterized by GC3 content. In human, the codons with poor-GC3 content will reduce mRNA stability (Hia et al. 2019). We found that the third bases of the 14 codons with higher usage frequency in SARS-CoV-2 were all A/T (Fig. 3F,G). In contrast, the 11 codons with higher usage frequently in the up-regulated genes ended with five G/C and six A/T (Fig. 3F,H). These results indicated that SARS-CoV-2 replication used the codons with poor-GC3 content, whereas the up-regulated genes in COVID-19 patients had no obvious preference for the third base of the codon. Previous studies have showed that the SARS-CoV-2 genome has a low GC content, and the codons ending in A or T were preferentially used (Kandeel et al. 2020). In agreement with this, we found that, compared with the human genome, the SARS-CoV-2 genome used the codons with poor-GC3 more frequently (Supplemental Fig. S3D). These results suggested that the increased expression of the AT3 codons in COVID-19 patients could reduce the stability of the host mRNAs and affect the synthesis of host proteins. Through time-course analysis, we found that the 14 codons related to SARS-CoV-2 had the highest abundance in the early stage of treatment, and their abundance gradually decreased as the treatment progressed (Supplemental Fig. S3E), indicating that, similar to the immune-responsive genes, the tRNA pool also returned to a normal state during therapy.

To investigate whether these GC3-poor codons were specifically up-regulated during SARS-CoV-2 infection, we compared the differences in codon usage and abundance under infection between SARS-CoV-2 and influenza A virus. FLUAV is also an RNA virus whose genome has a low GC content and prefers AT3 codons (Wong et al. 2010; Ahn and Son 2012; Pavon-Eternod et al. 2013), and a previous study indicated that, in human cells under FLUAV infection, 10 codons were up-regulated (Pavon-Eternod et al. 2013). We found that, among the 14 codons with both higher expression levels in COVID-19 patients and higher usage frequency in SARS-CoV-2 genome, six were also up-regulated in FLUAV infection, and their codon usage frequencies were also higher than the human genome (Supplemental Fig. S3F). Meanwhile, for the remaining eight codons whose abundance was not up-regulated in influenza A virus infection, seven of them were used more frequently in the SARS-CoV-2 genome (Supplemental Fig. S3F). These results indicated that, although the influenza A virus genome also had a preference for GC3-poor codons, eight codons were still specifically up-regulated during SARS-CoV-2 infection.

Microbial profiles of COVID-19 patients revealed by plasma cfRNAs

Some pathogenic microorganisms, like human mastadenovirus C (Saban-Ruiz and Ly-Pen 2020; Mehta et al. 2021), Haemophilus influenzae (Lansbury et al. 2020; Zhong et al. 2021), Candida glabrata (Arastehfar et al. 2020; Huang et al. 2020b; Posteraro et al. 2020), Klebsiella pneumoniae (Huang et al. 2020b; Lansbury et al. 2020; Arcari et al. 2021), and Pseudomonas aeruginosa (Huang et al. 2020b; Lansbury et al. 2020; Garcia-Vidal et al. 2021), have been detected in COVID-19 patients by sequencing and bacterial culture. According to the analysis of the cfRNA landscape, there were still 20% of reads that failed to map to the human genome (Supplemental Fig. S1A). We collected the unmapped reads and calculated the microorganism abundance in each sample using Kraken2 (Wood et al. 2019). We found that previously reported human mastadenovirus C (Saban-Ruiz and Ly-Pen 2020; Mehta et al. 2021), Klebsiella pneumoniae (Huang et al. 2020b; Lansbury et al. 2020; Arcari et al. 2021), and Candida glabrata (Arastehfar et al. 2020; Huang et al. 2020b; Posteraro et al. 2020) have higher abundance in COVID-19 patients compared with healthy donors (Fig. 4A). To further validate this result, we performed transcript assembly and visualized the assembled genes in the Integrative Genomics Viewer (IGV) (Fig. 4B–D; Robinson et al. 2011). In addition, the presence of the pathogenic microorganisms was also confirmed by clinical cultures and ultra-deep metatranscriptomic sequencing of matched respiratory specimens, including sputum or nasal secretions, based on our previous research (Zhong et al. 2021). Compared with healthy donors, the read coverage and the number of assembled genes were significantly increased, indicating that the microorganisms identified were bona fide rather than noise. One recent study found that there may be latent infections including Toxoplasma gondii infection in COVID-19 patients (Roe 2021). Consistent with the study, we also identified Toxoplasma gondii (Fig. 4A,E). These codetected microorganisms, especially pathogenic bacteria under SARS-CoV-2 infection, can aggravate the lung injury in COVID-19 patients, and microbial profiles of COVID-19 patients revealed by plasma cfRNAs open a novel angle in the identification of microorganisms and evaluation of disease status.

Figure 4.

Microbial detection in COVID-19 patients. (A) Box plot showing the expression levels (RPMs) of human mastadenovirus C, Klebsiella pneumoniae, Candida glabrata, and Toxoplasma gondii in healthy donors and COVID-19 patients. Each point represents the mean value of the RPM of a patient/healthy volunteer at all time points. (B–E) The IGV view showing cfRNA reads mapped to human mastadenovirus C genome (B), Klebsiella pneumoniae genome (C), Candida glabrata genome (D), and Toxoplasma gondii genome (E) in healthy donors and COVID-19 patients. (COV) All COVID-19 patients (including mild and severe COVID-19 patients), (control) all healthy donors. Asterisks indicate statistically significant differences; (***) P < 0.001.

Discussion

The rapidly spreading SARS-CoV-2 virus has had a great impact on human lives. At present, the mechanism of COVID-19 pathogenesis from the point of cfRNAs remains unclear. Here, we collected and sequenced plasma cfRNAs from both COVID-19 patients and healthy donors and obtained the cfRNA profiles of various RNA molecules. Among them, we identified 380 cfRNA molecules with increased expression levels in COVID-19 patients, and their expression levels were able to separate COVID-19 patients from healthy donors, including seven genes that could serve as potential COVID-19 biomarkers with AUC > 0.85. We observed enhanced antiviral (NFKB1A, IFITM3, and IFI27) and inflammatory responses in COVID-19 patients, especially neutrophil activation–related processes in severe patients, which was consistent with clinical data and previous studies (Huang et al. 2020a; Thierry and Roch 2020). Nevertheless, IFI27 was also able to distinguish influenza and noninfluenza flu-like illnesses (Tang et al. 2017) and was up-regulated during influenza virus, respiratory syncytial virus (Ioannidis et al. 2012), and human rhinovirus (Zhai et al. 2015) infections according to transcriptional analysis. IFITM3 can inhibit virus-triggered induction of type I interferon and was up-regulated to prevent influenza infection (Jiang et al. 2018). The specificity of some potential biomarkers in differentiating between SARS-CoV-2 infection and other respiratory viruses remains unstudied. Together with the time-course analysis, these results indicated that inflammatory responses were elevated in COVID-19 patients.

IL6/IL6R-mediated cytokine storm is a major feature of COVID-19 patients (Huang et al. 2020a). The protein level of IL6 was increased, whereas its mRNA level was not altered, indicating that up-regulation occurred at the translational level. Consistently, we indeed observed a global down-regulation of miRNA genes in both mild and severe COVID-19 patients. Further analysis revealed that the down-regulation of let-7 family members hsa-miR-451a, hsa-miR-23a-3p, and hsa-miR-23b-3p, which can directly target IL6/IL6R mRNAs (Iliopoulos et al. 2009; Aghaee-Bakhtiari et al. 2015; Liu et al. 2016; Huang et al. 2017; Sui et al. 2019), may promote the expression of IL6/IL6R in COVID-19 patients at the protein level. Meanwhile, an up-regulated lncRNA, GJA9-MyCBP, may function as a miRNA sponge to compete with IL6R/IL6 for the let-7 family. Therefore, the derepression of IL6/IL6R translation mediated by these noncoding RNAs may thereby enhance the cytokine storm in COVID-19 patients.

tRNAs are required for protein synthesis, and codon abundance determines translational efficiency and mRNA stabilities (Hia et al. 2019). The SARS-CoV-2 genome has a usage preference for AT3 codons, whereas the human genome prefers GC3 codons. We found that the tRNA pool was significantly changed, leading to the accumulation of SARS-CoV-2-biased codons and may thus facilitate SARS-CoV-2 replication. Meanwhile, this can also reduce the stability and translational efficiency of host mRNAs (Hia et al. 2019) and affect the synthesis of host proteins (Dilucca et al. 2020). In the future, the reduction of the expression levels of SARS-CoV-2-preferred tRNAs in host cells may be able to inhibit its replication due to insufficient tRNAs for its protein synthesis. On the other hand, codon deoptimization, a new technology that can attenuate the virus virulence by introducing the least-preferred codons (Nogales et al. 2014; Li et al. 2018), could be used as a novel strategy for SARS-CoV-2 vaccine development. Therefore, codon preference can also be considered in the design of future vaccines and drugs, which will help the development and design of vaccines to protect against COVID-19 (Ko et al. 2005; Lopes et al. 2017; Latanova et al. 2018).

COVID-19 patients had co-infected microorganisms, including Klebsiella pneumoniae (Huang et al. 2020b; Lansbury et al. 2020). In our data, we identified several reported and novel codetected microorganisms in COVID-19 patients. These codetected microorganisms, especially pathogenic bacteria under SARS-CoV-2 infection, may aggravate the lung injury and lead to an exacerbated inflammatory response, according to previous studies. Meanwhile, distinct signatures of microbial profiles among COVID-19 patients can be related to disease severity, according to our previous findings (Zhong et al. 2021). These results indicate that pathogenic microbial information can be revealed from the plasma of COVID-19 patients, providing additional information for COVID-19 therapy and monitoring, such as antibiotics selection and COVID-19 tracking (Lai et al. 2020).

As with the majority of studies, the design of the current study is subject to limitations. First, the sample size enrolled in this study is insufficient; people from other racial groups infected with different pathogens should be enrolled to confirm the findings in the future. Second, the relationship between cfRNA profile and clinical drug use was uncertain, which may influence disease process surveillance. Third, although the current study provides useful biomarkers discovered in the cohort of COVID-19 patients, the specificity of some potential biomarkers in differentiating between SARS-CoV-2 infection and other respiratory viruses is still lacking. Follow-up studies need to be conducted to explore the impact of treatment on cfRNA profile and the cfRNA characterization among different viral infections.

As a noninvasive approach, plasma cfRNA provides a good resource for COVID-19 research (Fig. 5), including potential biomarkers of disease severity (Fig. 5A), coding/noncoding gene regulatory mechanisms (Fig. 5B,C), and microbial community composition (Fig. 5D). Currently, multiple RNA sequencing methods have been applied to cfRNA analysis, including PALM-seq, Phospho-RNA-seq (Giraldez et al. 2019), small RNA-seq (Murillo et al. 2019), etc. Each platform has its unique advantages and disadvantages. Here, we adopt PALM-seq, which has been compared with multiple public data produced by SMARTer Seq, ScriptSeq, and small RNA-seq in our previous report (Yang et al. 2019). Compared with the traditional RNA library construction that requires a lot of blood, it is difficult to prepare large-scale cfRNA libraries with the high cost for large-scale cfRNA library preparation and low mapping rate. In contrast, the PALM-seq database construction method employs terminal modification and adds a 3′ polyadenylic acid and 5′ linker; thus, mRNA and noncoding RNAs (including long noncoding RNA, microRNA, tRNA, piRNA, and other RNAs) can be obtained in one library. These RNAs can be sequenced at a relatively low depth. In addition, PALM-seq enables the sequencing of the entire transcriptome in low-concentration RNA samples. This method allows us to perform more comprehensive cfRNA large-scale population research and clinical applications. Taken together, we reported the dynamic and unique cfRNA landscape of COVID-19 patients, which expands our knowledge of SARS-CoV-2 pathogenesis from a new perspective. This study also provides novel potential targets for the design of drugs or vaccines.

Figure 5.

Plasma cfRNA landscape provides resourceful solutions for COVID-19 research. (A) Activated anti-virus responses in both mild and severe patients may attenuate lung infection in COVID-19 patients. (B) Enhanced cytokine storm mediated by differentially expressed miRNAs and lncRNAs may aggravate lung injury in COVID-19 patients. (C) tRNA pool perturbation may facilitate SARS-CoV-2 infection which could be considered as a novel target for SARS-CoV-2 vaccine development. (D) Pathogenic microorganisms can be codetected in COVID-19 patients.

Methods

Patients and data collection

A total of 37 COVID-19 patients and eight age-matched healthy subjects were recruited from four Guangdong local hospitals from January 27 to February 26, 2020. Whole blood was collected using BD Vacutainer spray-coated K2EDTA tubes in the morning. No participants were fasting and the time interval between the patient's blood collection and treatment was at least 6 h. In addition, we collected clinical data including blood routine, blood biochemistry, blood coagulation function, interleukin 6, and days postadmission (Supplemental Table S1). After completing routine blood tests in the clinical laboratory, the leftover surplus blood from participants who signed informed consents was used for subsequent cfRNA study. Those samples with sufficient volume (≥1 mL/sample) were retained and transferred to a Biosafety Level 2 (BSL-2) laboratory at ambient temperature. Plasma separation was performed within 3 h after blood collection by using the two-step separation procedure (Chiu et al. 2001): 1600g for 10 min at 4°C and 16,000g for 10 min at 4°C. Opaque or turbid plasma at the first or second centrifugation step was excluded to prevent the rupture of cellular DNA. Then, RNA was extracted immediately using the HiPure Liquid RNA Mini kit (Magen) at this BSL-2 laboratory. A NanoDrop One instrument was used to check RNA concentration and quality. Parallel treatment and analysis of plasma from healthy donors were used for quality control and comparative analysis. The extracted RNA was stored at −80°C until used. The RNA library was constructed using the PALM-seq (Yang et al. 2019) protocol which can capture plasma mRNA, lncRNA, and small RNA with high complexity. The steps are as follows: (1) 3′ polyadenylation was added by Escherichia coli poly(A) Polymerase (PAP); (2) 5′ adaptor was ligated by T4 RNA Ligase 1; (3) rRNAs or other RNAs of no interest were consumed and processed by RNase H and DNase I; (4) RNA with 3′ polyadenylation and 5′ adaptor was reverse-transcribed by oligo(dT) with 3′ adaptor and amplified by PCR. We quantified the concentration of cfRNA library products using a Qubit 3.0 fluorimeter (Supplemental Table S1). After quality control, the high-quality library products were sequenced on the DNBSEQ platform (BGI) in single-end 100 bp.

CfRNA data analysis

The fastp (version 0.20.1) (Chen et al. 2018) tool was used to check the base quality, connection condition, GC content, sequence length, and repeat sequence, and to remove poly(A) tails, sequencing adapters, and short reads (<17 bp) with parameters “qualified_ quality_phred = 5, n_base_limit = 15, unqualified_percent_limit = 50”. After removing poly(A) tails, sequencing adapters, low-quality bases and human ribosomal RNA, clean sequencing reads were mapped to the human reference genome (GRCh38.p12) and the SARS-CoV-2 genome (GCF_009858895.2) using HISAT2 (version 2.0.4) (Kim et al. 2015) with parameters “-k 1 -q ‐‐novel-splicesiteoutfile ‐‐no-unal ‐‐dta ‐‐un-gz ‐‐rna-strandness F”. Then, the expression levels of each gene were calculated by the transcripts per kilobase of exon model per million mapped reads (TPM) using StringTie (Pertea et al. 2016) with parameters “-t -C -e -B -A ‐‐fr” based on the result of HISAT2.

Differential gene expression analysis

A 1.5-fold variance in expression levels and a P-value <0.05 were used as cutoffs to define differentially expressed genes. The P-value was calculated using R (R Core Team 2021) software (DESeq2) (Love et al. 2014).

Hierarchical clustering analysis

Hierarchical clustering (HC) was performed with up-regulated and down-regulated genes in severe and/or mild COVID-19 patients, compared with healthy donors. The distance used was Euclidian, and the cluster method was Ward's method (ward.D2).

Gene Ontology analysis

Gene Ontology analysis was performed using the R (R Core Team 2021) software (clusterProfiler) (Yu et al. 2012). The annotation Dbi R package org.Hs.eg.db was used to convert gene symbol to gene id. A P-value <0.05 was taken to indicate statistical significance.

Time-course analysis

In total, gene expression data from 16 severe and 14 mild patients were used for the time-course analysis (Supplemental Table S1). First, we identified 682 and 431 genes differentially expressed between three stages in severe and mild patients, respectively. Then, we used the R (R Core Team 2021) package (mfuzz) (Futschik and Carlisle 2005) to perform time-course analysis with the fuzzification parameter calculated by the function mfuzz::mestimate.

miRNA target gene prediction

The R (R Core Team 2021) package (miRNAtap) was used to predict the target genes of 41 down-regulation miRNAs based on the annotation package of human genes org.Hs.eg.db. We used the getPredictedTargets function to predict the mRNA targeted by the miRNA; the parameters are species = “hsa” and method = “geom”. Experimentally verified miRNAs and mRNAs regulatory information were obtained from the miRTarBase database (Huang et al. 2020c) (http://miRTarBase.cuhk.edu.cn/).

Codon abundance analysis

The anticodon sequences of tRNAs were first converted to the corresponding codon sequence, and the expression levels of tRNAs with the same codon were merged. The P-value was calculated using R (R Core Team 2021) software (DESeq2) (Love et al. 2014). A twofold variance in expression levels and a P-value <0.05 were used as cutoffs to define differentially expressed codons.

Nucleotide composition analysis

The nucleotide composition of the codons was calculated by WebLogo 3 (http://weblogo.threeplusone.com/create.cgi) (Crooks et al. 2004).

Codon usage frequency analysis

The codon usage frequency of human, SARS-CoV-2, and FLUAV was calculated by the number of each codon in the coding sequence divided by the total codon number.

Receiver operating characteristic analysis

The ROC curve of 380 up-regulated genes and 138 up-regulated tRNAs in COVID-19 patients was calculated using R (R Core Team 2021) software (ROCR). The performance function was used to calculate the area under the curve (AUC), and the parameter is measure = “auc”. First, the sample was divided into two groups. One group was the COVID-19 (COV) group, including severe and mild patients, and the other group was a control group composed of healthy donors. Then, the R package ROCR was used to calculate the AUC values for each up-regulated gene and 138 up-regulated tRNAs in COVID-19 patients. The AUC score of a gene or 138 up-regulated tRNAs showed the probability of whether the gene or tRNAs will distinguish the two defined sets of data. The higher the AUC, the better the effectiveness of distinguishing the two sets of data. Finally, we used 0.85 as the threshold to screen out genes with AUC > 0.85, and these genes were the predicted biomarkers for COVID-19.

Logistic regression analysis

The SelectFromModel class and LogisticRegression class in the sklearn module were used for feature selection. The cross_val_score function was used for cross-validation. The samples were divided into two groups. One group is the COVID-19 (COV) group, including severe and mild patients, and the other group is a control group composed of healthy donors.

Identification of codetected microorganisms

PALM-seq (Yang et al. 2019) reads that failed to map to the human genome were collected and mapped to the microbial database using Kraken2 (version 2.0.8-beta) (Wood et al. 2019) with parameters “‐‐threads 4 ‐‐gzip-compressed”. The abundance of each microorganism was calculated by reads per million mapped reads (RPM). The P-value was calculated using the Wilcoxon test, and the P-value <0.05 was taken to indicate statistical significance. The fold changes of the microorganism abundance in mild and severe patients compared with healthy donors were calculated respectively. A twofold variance in abundance and a P-value <0.05 were used as cutoffs to define differentially expressed microorganisms.

Gene assembly of codetected microorganisms

The reads that failed to map to the human genome were collected and then mapped to the human mastadenovirus C genome (GCA_006433735.1), Klebsiella pneumoniae genome (GCF_000240185.1), Candida glabrata genome (GCF_000002545.3), and Toxoplasma gondii genome (GCF_000006565.2) using HISAT2 (version 2.0.4) (Kim et al. 2015) with parameters “-k -p -x ‐‐no-spliced-alignment -q ‐‐no -unal ‐‐dta”. Transcript assembly was performed using StringTie (Pertea et al. 2016) with parameters “-p ‐‐rf -m -t”, and then the assembly results of different samples were combined and assembled according to control and disease using Cuffcompare with parameters “-r -R -s -i”.

Data access

All raw and processed sequencing data generated in this study have been submitted to the CNGB Sequence Archive (CNSA; https://db.cngb.org/cnsa/) (Guo et al. 2020) of the China National GeneBank DataBase (CNGBdb) (Chen et al. 2020b) under accession number CNP0001306.