Biomarkers and Immune Repertoire Metrics Identified by Peripheral Blood Transcriptomic Sequencing Reveal the Pathogenesis of COVID-19

The work was performed in line with the Declaration of Helsinki principles for ethical study. The study was approved by the Medical Ethics Committee of the Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China [Clinical Ethical Approval No. (2020) IEC-A251].

The RNA-seq data of acute and healthy control were retrieved from the NCBI-GEO dataset, with ID PRJNA639275 (27). RNA-seq raw data were cleaned using Fastp (v0.20.1) and then mapped to human genome (GRCh38) using HISAT2 (v2.2.0). FeatureCounts tool of Subread package (v2.0.0) was used for counting the number of reads aligned to genes, with the help of GENECODE gene annotation (v34). Then analysis of differentially expressed genes was performed with DESeq2 R package (v1.30.0). Genes were categorized as differentially expressed with the criteria: adjusted p-value (pAdj) <0.05 and |log2FoldChange| >1 for upregulated genes and the adjusted p-value <0.05 and |log2FoldChange| <−1. Gene counts were normalized by DESeq2 and converted to log₂ (normalized counts + 1) format. Then batch effect was removed by limma R package (v3.46.0) and the result was taken as the normalized gene expression level.

GO and KEGG enrichment analyses of differentially expressed genes (DEGs) were performed using clusterProfiler R package (3.18.0). The Benjamini and Hochberg method was used to adjusted p-value. Significantly enriched GO term and KEGG pathway were selected under the criteria: adjusted p-value <0.05.

We applied MiXCR (v.3.0.12) to RNA-seq data from peripheral blood mononuclear cells (PBMCs) for reconstructing B-cell receptor (BCR) and T-cell receptor (TCR) reads. To minimize batch effects and address the unevenness of sequencing depth, sequencing data volume (the number of reads multiplied by the length of reads) was kept at the same level. The targeted sequencing of the TCR beta chain CDR3 to distinguish T-cell clones, frequency, and antigenicity is feasible in TCR repertoire identification (28, 29).

BCR and TCR clonotypes were extracted using MiXCR (v.3.0.12). "Analyze shotgun" command was used with parameters "–impute-germline-on-export" and "–only-productive". The post-analysis of BCR/TCR repertoire was conducted using VDJtools (v1.2.1).

The number of unique clonotypes per sample was used to indicate richness. Unique clonotype was strictly defined by exact match for receptor nucleotide sequence (CDR3 sequence, V/J segments, and hypermutations). Shannon entropy was used to represent diversity. Shannon entropy is expressed as:

Where N represents the total number of clonotypes in sample X and p_xi refers to the frequency of clonotype i. Generally speaking, more unique clonotypes and a more even distribution of clonotypes will increase the diversity of a given repertoire.

1-Pielou' evenness was used to represent clonality which ranges from 0 to 1, and zero clonality means that each clonotype has the same frequency. Higher clonality often means that there is a group of highly expanded clonotypes with dominant frequencies.

Pielou' evenness is expressed as:

Where H represents Shannon entropy and H_max is calculated as ln(N). (N represents the total number of clonotypes in a sample).

The clonality is expressed as:

In CDR3 length analysis, the samples in each group were analyzed all together. Length distribution was described by cumulative distribution function (CDF).

We used a logistic regression model as the classifier. The logistic regression function is expressed as follows:

Where f₁ to f_n represent different features. b₀ indicates bias and b₁ to b_n denote the weight of these features. A sigmoid function was applied to transform logit into a score ranging between 0 and 1. The sigmoid function is expressed as follows:

As a Python module integrating classical machine learning algorithms, Sklearn (v0.23.1) was used to implement our model. To evaluate the performance achieved by each combination of features, a three-fold cross-validation was performed. The samples were randomly split into three sets and each set was called a fold. For each fold, it was taken as a test set and the remaining two folds were treated as a training set. Then, this progress was repeated until all folds had served as the test set. The performance of a model was characterized by the mean accuracy and mean AUC of results tested on each fold. For each combination, a ROC curve was outputted according to the result of the prediction. The average AUC of the curves was compared and the final output combination was determined with the highest AUC.

SARS-CoV2 S1 protein (Sino Biological Inc., BJ) was coated on a high-binding 96-well plate (Thermo Scientific) at 4°C overnight. Plates were blocked with 5% non-fat milk in phosphate-buffered saline (PBS) for 1 h at room temperature, followed by incubation with 1:100 diluted plasmas in dilution buffer (PBS, 2% non-fat milk, and 0.05% Tween-20). A 1:4,000 dilution of horseradish peroxidase (HRP)-conjugated mouse anti-human IgG, IgM, and IgA antibodies (BaiaoTong Experiment Center, LY) was added and incubated for 1 h at room temperature. Wells were washed six times between each step with 0.05% Tween-20 in PBS (PBST). Finally, wells were developed using tetramethylbenzidine substrate (Beyotime Inc., WH) and were stopped by the addition of stop solution, followed by reading at 450 and 570 nm. The sample OD was calculated as OD450 subtracted by OD570.

Statistical analyses were performed with GraphPad Prism software (Version 6.0, GraphPad Software Inc.). Throughout the study, n refers to the number of subjects where every subject is one data point. Unpaired two-group comparisons were done with the Mann–Whitney U test. In figures, ****P < 0.0001, ***P < 0.001, **P < 0.01, and *P < 0.05.

A total of 62 COVID-19 convalescents, 16 acute patients, and 17 healthy donors were involved in this study. The information of subjects including age, sex, disease severity, and time from symptom onset or a negative RT-PCR test to sample collection is summarized in Table 1. Of note, the samples of healthy controls were collected prior to the outbreak of COVID-19. Acute patients with COVID-19 were all hospitalized, while 60 of the 62 convalescents were hospitalized. The percentages of hospitalized patients in the acute infection group and convalescent group were 100% and 96.7%, respectively. The median time from symptom onset to sample collection for acute patients was 14 days, with an interquartile range of 9–16 days. In the convalescent group, the median time from a negative RT-PCR test to sample collection was 151.0 days, with an interquartile range of 135–162 days. The majority of these convalescents (51 in 62; 82.2%) were moderate cases, and only 2 were mild cases, and 9 were severe cases (Table 1).

Ficoll density gradient centrifugation was performed to isolate the PBMCs, followed by the construction of RNA library and high-throughput sequencing (Figure 1A). Quality control was applied for the RNA-seq data, and batch effect was minimized by including a pooled mixture of all samples as standard control.

In order to investigate the peripheral alteration during and after SARS-CoV-2 infection, we compared the transcriptomic differences in PBMC samples from acute COVID-19 patients, convalescents, and healthy donors through DEG analysis. The combined DEG list is shown in the. 2

The heat map of the DEG profiles was built (Figure 1B), through which the profile of acute COVID-19 patients was found to be clearly different from that of convalescents and health donors (HDs). Of note, the DEG profile of convalescents was comparable to that of HD, suggesting that the overall gene expression pattern of acute patients was significantly altered compared with that of HDs, and almost returned to the normal level of healthy individuals in the convalescent phase.

To further study the transcriptomic changes induced by SARS-CoV-2 infection regarding the specific cellular biological processes, we performed a volcano plot (Figure 1C) followed by gene functional enrichment analysis (Figure S1) to identify differentially expressed genes and enriched pathways.

In the comparison between acute COVID-19 patients and HDs, 668 genes were upregulated and 230 genes were downregulated in acute COVID-19 patients. Whereas in the comparison between convalescents and acute patients, 475 DEGs were downregulated in convalescents, compared with 91 genes that were upregulated. When the DEGs of convalescents and HDs were compared, sparse DEGs were noted and the fold change of gene expression was not evident, indicating the resemblance between these two cohorts from the transcriptomic perspective.

Interestingly, the upregulated genes in acute patients compared with those in HDs were enriched in multiple biological processes, which were also found to be downregulated in convalescents compared with those in acute patients, as revealed by GO term and KEGG pathway enrichment analysis (). Of note, the viral infection-induced alterations in PBMCs of convalescents were enriched in gene modules including "complement activation," "humoral immune responses," "phagocytosis," and "adaptive immunity based on somatic recombination." In addition, the upregulation of a series of B-cell-related gene modules was also observed, including "immunoglobulin receptor binding" and "antigen binding." These findings suggested that the expression of genes related to humoral immunity may be upregulated in the acute phase of SARS-CoV-2 infection, thereafter declining and finally returning to the ground level in convalescents. To conclude, our data provided a quantitative and global insight into the dynamic changes of immune-related genes, such as those implicated in B-cell response during SARS-CoV-2 infection, especially in the convalescent phase. 3

To study the alteration of the B-cell and T-cell repertoire during and after SARS-CoV-2 infection, we retrieved the complementarity determining region 3 (CDR3) variable region of the beta chain of the TCR and heavy chain of the BCR from the RNA-seq library. B-cell diversity, as estimated by unique CDR3 usage, was significantly higher in acute COVID-19 patients compared with that in convalescents (P < 0.0001) and HDs (P < 0.0001) (Figure 2A). Next, richness, a measurement used to quantify the number of unique V(D)J rearrangements, was evaluated and also found to be significantly higher in acute COVID-19 patients than in convalescents and HDs. Finally, clonality, an indicator used to characterize clonal amplification and represents the evenness of the repertoire, was higher in acute patients than in healthy controls and recovered people. To better explain the meaning for quantification of B-cell diversity, richness, and clonality, a toy model was used and displayed in Figure S2. In sum, the highest clonality, diversity, and richness of B-cell repertoire were observed in acute COVID-19 patients, highlighting a biased and uneven B-cell repertoire developed in peripheral blood during the acute phase of SARS-CoV-2 infection, which may reveal the dynamic nature of B-cell response during SARS-CoV-2 infection.

The overall T-cell metrics for these subjects were also calculated according to the CDR3 sequence of the beta chain of the TCRs (TRB). Notably, the highest level of TRB diversity was observed in HDs (Figure 2B), whereas the TRB diversity of acute COVID-19 patients and convalescents was significantly lower than that of HDs (P < 0.0001, t-test). However, the acute patients and convalescents had significant higher T-cell clonality than HDs, indicating the clonal expansion of antigen-specific T cells and the shift of T-cell repertoire. The richness of TRB sequence, which represents the distribution of unique clonotypes, showed a similar trend as the diversity. Together, lower diversity and evenness as well as higher clonality of T cells in acute and convalescent patients may indicate the durability of T-cell response to SARS-CoV-2 infection.

We next sought to compare the profiles of TCR gene usage between the three groups. In this regard, significantly distinct profiles were observed in TRB genes (Figure 2C), but were less pronounced in TRA (alpha chain of the TCRs) genes (Figure S3). An overrepresentation of TRBJ2-7 was identified in convalescents, whereas an overrepresentation of TRBJ1-1 and TRBJ2-1 was observed in acute patients. Interestingly, we also discovered TRBV28/TRBJ2-7 as the top V-J pair in convalescents. The distinct usage of V(D)J genes between acute patients and convalescents could be explained by a repertoire shift reflecting the transition from acute T-cell response to memory T-cell response.

Previous studies had shown that TCR structural difference, such as CDR3 length, might contribute to the differentiation of T cells (30). To test this issue, we analyzed the CDR3 length of BCR and TCR in the three cohorts and found that B cells from acute patients tended to have a longer IGH (heavy chain of immunoglobulin) CDR3, whereas no differences were observed in T cells (Figure 2D and Figure S4), suggesting that heterogeneous BCR recombination might be induced mainly in IGH upon SARS-CoV-2 infection.

Taken together, these results highlighted the complex and imbalanced nature of host immune responses to SARS-CoV-2 infection, in which B-cell response may be substantially induced in the acute phase and return to the ground level in the convalescent phase, whereas T-cell response could persist for a long time even without the presence of specific antigen in the convalescent phase.

Next, a computational pipeline was developed based on the transcriptomic data of our cohort to identify the potential biomarker or biomarker combination that could accurately distinguish between different disease courses of COVID-19. More importantly, we aimed to distinguish the convalescents from the unexposed HDs, since some convalescents who were asymptomatic could not be successfully identified by nucleic acid testing (NAT) or antibody testing.

Firstly, all transcriptomic data were normalized on indexes including repertoire metrics and DEGs before machine learning algorithm was applied. Twenty-nine biomarkers consisting of 20 repertoire metrics and 9 highly ranked DEGs chosen from 1,241 DEGs were taken as candidate biomarkers to generate about 146,595 groups of candidate biomarker combinations. Samples were randomly divided into two groups, with 70% and 30% of the samples allocated to the training set and the testing set, respectively. Cross-validation was preformed using these datasets to determine the preferred biomarker combination with the highest area under the curve (AUC) value (31 –33). For the model training in this step, multivariate logistic regression (LR), a widely accepted machine learning algorithm, was performed (Figure 3A). This approach resulted in the identification of five independent biomarkers which could be used to predict disease outcome and distinguish between acute patients, convalescents, and HDs (Figure S3). These biomarkers include the differentially expressed gene C1D nuclear receptor corepressor (C1D), Alpha-2-adrenergic receptors (ADRA2A), spindlin family member 3 (SPIN3), ubiquitin C-terminal hydrolase L1 (UCHL1), and Glutathione S-Transferase Mu 3 (GSTM3). To determine the reliability of these biomarkers, confusion matrices (Figure 3C) were built and showed that various diseases courses of these samples could be distinguished with ideal accuracy with these biomarkers. Moreover, a best biomarker combination was also determined, including B-cell diversity, T-cell diversity, and clonality, as well as the expression levels of UCHL1 and C1D. This combination included both DEGs and repertoire features and thus resulted in optimal accuracy of the model. The AUC of acute, convalescent, and HD reached 0.99 in this model (Figure 3B). Moreover, the principal component analysis (PCA) (Figure 3D) demonstrated that the different clusters of samples were clearly segregated by this combination, further suggesting the reliability of this prediction model.

Thus, our findings demonstrated that the clinical outcome of a patient can be predicted by a combination of repertoire signatures and DEGs through machine learning approaches, thereby providing insight into the classification and prediction of COVID-19 development and invoking timely and precise clinical advice for the treatment of COVID-19.

Among the aforementioned biomarkers, UCHL1 is particularly interesting given the following reasons: i) it is a biomarker frequently used for indicating brain damage induced by viral infections; ii) the encoded protein can be tested in the blood sample, making it clinically applicable; iii) it is one of the five "predictors" identified by our machine learning approach. Therefore, we next focused on this molecule. Firstly, we confirmed the normalized gene expression profiles of UCHL1 measured by RNA-seq in each cohort with or without COVID-19 (Figure 3E). After normalization, the expression level of UCHL1 was still found to be significantly higher in acute patients than in convalescents and HDs.

Therefore, we next sought to explore the association of clinical outcomes with UCHL1 levels at the time of sampling (Figure 3F). For patients with severe symptoms or admitted to the ICU, the UCHL1 level was negatively associated with the time from symptom onset to sample collection (Corr = −0.554, P = 0.0049). A tendency of negative correlation was also noted when all acute patients were included for analysis; however, this correlation failed to reach statistical significance (Corr = −0.173, P = 0.344).

Together, we found that UCHL1 expression used as a biomarker in the peripheral blood can predict the disease development of COVID-19, although the potential mechanism by which this molecule participates in the pathogenesis of COVID-19 requires further investigations.

To further elaborate the role of UCHL1 in COVID-19, we reviewed human studies investigating the use of serum UCHL1 as biomarker for brain injury. In a multicenter cohort study including 206 patients and 175 healthy controls, UCHL1 level was effective in discriminating between TBI (traumatic brain injuries) patients with and without intracranial lesions on a CT scan and was correlated with functional 3-month outcomes (34). In another study with 273 mild and moderate TBI patients, similar results were observed and the combination of GFAP and UCHL1 was sensitive for identifying a positive CT scan (35). Another observation cohort study on HIE (hypoxic–ischemic encephalopathy) patients (neonate patients = 16, controls = 11) also found that UCHL1 was elevated in HIE neonates and associated with cortical injury (36). As a protein biomarkers for central nervous system damage, UCHL1 has been detected in diverse brain injuries including ischemic/hemorrhagic stroke (37, 38), TBI (39, 40), Parkinson's disease (41), cardiac arrest (42), and seizures (43).

The role of human coronaviruses in nervous system damage has been underestimated. Nevertheless, recent evidence suggests that the development of SARS-CoV-2 infection could induce neurological disorders involving both the central and peripheral nervous systems in COVID-19 patients worldwide (44 –48); however, the pathogenesis and the sequelae of the damage are still poorly understood. Our transcriptional analysis suggests that the elevated serum UCHL1 protein levels could be a hallmark of severe illness in acute COVID-19 patients; theoretically, this protein might be involved in opportunistic direct attack on the brain or systemic secondary insult including thrombosis, hypoxemia, or autoimmune response. Therefore, future clinical investigations are warranted to build the connection of this molecule with severe illness of COVID-19 at the protein level given the accessibility of blood sample and the potential clinical relevance.

The appearance of T-cell diversity in the final biomarker combination which can accurately predict disease outcomes suggests that the structure of T-cell repertoire might contribute to anti-virus response during SARS-CoV-2 infection. To validate this speculation, the S protein-specific IgG levels of blood samples of convalescent individuals (n = 62) were measured by enzyme-linked immunosorbent assay (ELISA) (Figure 4A and Table S1), and correlative analyses were performed to determine their relationship with the immune repertoire metrics extracted from transcriptional analysis. These analyses revealed that the S protein-specific IgG levels were found to be positively correlated with T-cell clonality and diversity (r = 0.27; P = 0.0089 and r = −0.18; P = 0.0912, Spearman rank correlation) (Figure 4B). It is noteworthy to mention that even from patients who had been discharged for more than 200 days, the S protein-specific IgG was still detectable with a relatively high titer, accompanied by the persistence of memory T-cell response. Surprisingly, the S protein-specific IgG levels were not associated with those B-cell repertoire metrics. Nonetheless, the lack of association between antibody levels and B-cell response has been reported by other groups, which requires further investigation (discussed later).

We next sought to determine the kinetics of the T-cell proliferation. For this purpose, time of sample collection was incorporated into our analysis. The T-cell diversity was negatively correlated with days after recovery (r = −0.334; P = 0.029, Pearson correlation) (Figure 4C). However, the clonality of T cells exhibited a growing trend over time (r = 0.371; P = 0.014, Pearson correlation). The opposite trends of diversity and clonality could be explained by the clonal expansion of specific types of T cells, suggesting the prolonged maturation of T-cell responses against SARS-CoV-2 infection.

To conclude, the correlation between T-cell expansion and the SARS-CoV-2-specific antibody responses highlighted the contribution of T cells to the humoral immune responses. Besides, the persistent presence of T-cell response and virus-specific IgG antibody supports the notion that immune memory was maintained in convalescent COVID-19 patients for up to 6 months.