Analyzing master regulators and scRNA-seq of COVID-19 patients reveals an underlying anti-SARS-CoV-2 mechanism of ZNF proteins

Clinical data source

The clinical data comes from 3044 COVID-19 patients who were hospitalized in Wuhan Huoshenshan Hospital from 4 February 2020 to 13 April 2020. All COVID-19 patients were diagnosed according to the ‘New Coronavirus Pneumonia Diagnosis and Treatment Plan’ (7th edition) issued by the National Health Commission of China. This study has been approved by the Research Ethics Committee of Wuhan Huoshenshan Hospital, China (HSSLL011). Written informed consent was obtained from each patient.

Public data source

The transcriptome data of 430 SARS-CoV-2 positive nasopharyngeal swabs and 54 negative controls are obtained from the GEO database (GSE152075) [14]. The PBMC and bronchoalveolar lavage fluid (BALF) single-cell transcriptome data comes from GSE149689 [15] and GSE145926 [16].

Differential expression analysis and MR inference

The edgeR R package is used to calculate differentially expressed genes [17, 18]. ARACNe and VIPER algorithms are combined to identify MRs [19, 20]. ARACNe algorithm is a method based on mutual information theory, which can infer the interaction between transcription factors (TFs) and target genes [19]. VIPER was used to infer the enrichment of each regulatory protein regulator in certain gene expression signature (GES) composed of full transcriptome data of healthy and COVID-19 patients [20]. These two algorithms have proved to be very effective for identifying the MR in the host’s response to viral infection [21, 22]. A method called Leading-edge Gene Set Enrichment Analysis implemented in the VIPER package was performed to identify the core target genes of MRs driving the enrichment on GES [20].

Gene function analysis and gene set activity calculation

The clusterProfiler package was used for gene ontology and KEGG pathway enrichment analysis [23]. The single-sample gene set enrichment analysis algorithm (ssGSEA) in the GSVA package was used to calculate the overall activity of ZNF proteins [24]. This algorithm can evaluate the overall activity of a predetermined gene set in a single sample based on the entire transcriptome expression. The 456 ZNF protein coding genes matched by the grep function were used as the predetermined antiviral gene set of ZNF proteins.

Single cell data analysis

The Seurat (version 3.0) was used for single-cell transcriptome data quality control, filtering, standardization and subsequent analysis [25]. All inclusion criterion and parameters for cell quality control and dimension reduction clustering were consistent with the data literature source [15, 16]. The canonical correlation analysis and anchors algorithm were used to remove batch effect [25] and high-dimensional cell data were visualized through uniform manifold approximation and projection (UMAP) [26]. We used cell scores to assess how well individual cells express a certain set of predetermined genes via AddModuleScore function in Seurat [25, 27].

Flow cytometry assay of lymphocyte subsets

A total of 548 individuals from 3044 COVID-19 patients were enrolled in this study, and measured with T cell, B cell and NK cell using flow cytometry. A total of 125 patients of them were measured with 30 types of lymphocyte subsets. The peripheral blood was collected from patients and all samples were tested within 6 h after being obtained. Briefly, 100 μL of fresh whole blood was incubated in 2 mL of VersaLyse (Beckman Coulter Life Science) between 20 and 30°C for 15 min to lyse erythrocytes, then being washed with 3 mL of 1 × PBS. Thereafter, the cell pellet was resuspended in 500 μL of 1 × PBS containing 0.8% IOTest 3 Fixative Solution. The above samples were ready for acquisition. To measure T cell, B cell, NK cell and lymphocyte subsets, stained all the 10 single color tubes from a single pouch of the Compensation Kit provided in the IM DuraClone IM cell subsets Tube, 25 tests (Table S1, see Supplementary Data available online at http://bib.oxfordjournals.org/), RUO with venous blood, measured by multiple-color flow cytometry according to the manufacturer’s instructions (Backman Coulter Life Science). The data were evaluated using the Kaluza Software of Beckman.

The infection rate and prognosis of different genders

As shown in Figure 1, the infection rate of males is only 1.58 percentage points higher than females (Figure 1A), but males have a lower survival rate (Figure 1B) and a higher critical disease classification (Figure 1C) as well as a higher intensive care unit (ICU) admission ratio than females (Figure 1D). The death age for males is mainly at 60 ~ 90 years (Figure 1E), but 70 ~ 90 years for females (Figure 1F). The mortality rate of males in the same age group is about twice higher than females (Figure 1E and F), whereas only slight differences exist in the infection rate of different genders in the same age group except for the minimum or maximum age group with few subjects (Figure 1G and H).

Identifying MRs of different gender COVID-19 patients

We applied the ARACNe algorithm to, respectively, obtain 113 676 and 103 621 transcriptional interactions for male and female patients and found that 839 TFs regulate 9427 target genes in males and 864 TFs regulate 9565 target genes in females. We next used the VIPER algorithm to identify the top 50 proteins as MRs for males and females, respectively (Figure 2A and B). We found that males and females share 16 common MRs, whereas 34 MRs are different (Figure S1, see Supplementary Data available online at http://bib.oxfordjournals.org/).

GO annotation of MRs of different gender COVID-19 patients

As shown in Figure 3A, these MRs unique to male patients mainly participate in lymphocyte differentiation, regulation of leukocyte differentiation and T cell differentiation, but females are involved in steroid hormone-mediated signaling pathway and glucocorticoid receptor signaling pathway, whereas the 16 common MRs are enriched in defense response to virus, response to type I interferon, regulation of viral process, indicating that these common MRs involve in antiviral infection. Remarkably, the expression extent of 16 common MRs existed obvious differences between different gender patients (Figure 3B), and the fold change of STAT2, ZNFX1, ZNF121, ZNF347, ZNF562, AFF1 and SP1 of females is higher than males, respectively (Figure 3B), as well as the average protein activity of 16 common MRs in females is higher than males (Figure 3C).

Female patients express more ZNF proteins and higher activities than males

Previous studies indicated that STAT1/STAT2 could activate the expressions of many ZNF proteins (Figure S3, see Supplementary Data available online at http://bib.oxfordjournals.org/) [28, 29]. Here, we demonstrated that mRNA expression levels and protein activities of STAT1 and STAT2 after SARS-CoV-2 infection are significantly upregulated than controls, in particular the upregulated level of SATA2 in female patients is higher than males (log₂FC: 1.25 versus 2.14) (Figures 3B and S2, see Supplementary Data available online at http://bib.oxfordjournals.org/). We also found that female patients have 2 times more ZNF proteins than males among the top 50 MRs, and in addition to 5 shared ZNF proteins, females have 11 unique ZNF proteins, but only 3 in males (Figure 4A), as well as the overall protein activity of ZNF proteins is higher than males (Figure 4B). Moreover, 50 ZNF proteins were significantly up-regulated in females (log₂FC > 1, FDR < 0.05), but only 12 in males, while these down-regulated ZNF proteins in females were equal to males (15 versus 15) (Figure 4C and D). And that the activity score value of ZNF proteins in females was higher than males (Figure 4E).

Total 43 102 single-cell data were further clustered into 27 clusters (Figure S4, see Supplementary Data available online at http://bib.oxfordjournals.org/). We used typical genes as immune cell markers to identify 12 main cell populations (Figures 4F and G and S5, see Supplementary Data available online at http://bib.oxfordjournals.org/). We found that male and female patients share the same immune cell types (Figure 5A), but females have much higher proportions of CD4 T cells, CD8 T cells and fewer monocytes than males (Figure 5B). Especially, the overall ZNF protein activity of females is higher than males whether in mild or severe patients (Figure 5C and D). Among 12 cell subgroups, the ZNF protein activity of mono DCs, GC B cells, CD14 monocytes and NK cells in mild female patients is, respectively, significantly higher than mild males (Figure 5E), in particular the ZNF protein activity differences of RBC, NK cells, GC B cells, CD8 T cells, CD4 T cells and CD14 monocytes are more obvious between severe female and male patients (Figure 5F). These above results indicated that female patients could express more ZNF proteins and stronger transcriptional activities than males do.

The activity of ZNF protein is related to anti-SARS-CoV-2 infection

Interestingly, the overall ZNF protein activity of COVID-19 patients was significantly higher than controls, and the ZNF protein activity of patients with high viral load was higher (Figure 6A). Remarkably, the increased ZNF protein activity was significantly negative correlation with the viral load of patients (Figure 6B), in particular the ZNFX1 is more significantly negative correlation with viral load (Figure 6C).

We further detected the overall ZNF protein activity in single cell level and found that healthy people and COVID-19 patients have the same cell subpopulation distribution (Figure 6D), but the proportion of CD4 and CD8 T cells in the serum of severe patients is, respectively, lower than healthy individuals and mild patients, and the monocytes are higher than healthy or mild ones (Figure 6E). The order of overall ZNF protein activity is healthy people > mild patient > severe patient (Figure 6F). We also found that the ZNF protein antiviral activity of CD4 T cells, CD8 T cells, CD14 and CD16 monocytes in severe patients are lower than healthy people and mild patients (Figure 6G), whereas the anti-SARS-CoV-2 activity of the other 8 cell subgroups has no obvious differences (Figure S6, see Supplementary Data available online at http://bib.oxfordjournals.org/). Remarkably, the ZNF protein activity is positively correlated with the abundance of 19 immune cells, especially T cells (Figure 6H), and it also is positively correlated with the abundance of T cells in scRNA-seq (Figure S7, see Supplementary Data available online at http://bib.oxfordjournals.org/).

ZNF proteins as TFs participate in anti-SARS-CoV-2 role

Interestingly, male and female patients shared five same ZNF TFs (Figure 4A), but their target genes were not identical, and females had more target genes than males (Figures S8 and S9, see Supplementary Data available online at http://bib.oxfordjournals.org/). Go enrichment analysis demonstrated that these target genes of male and female patients can significantly enrich in ‘interferon response’, ‘antiviral response’, ‘regulation of viral replication and cycle’ (Figure 7A) as well as a wide range of ‘RNA binding ability’ molecular function (Figure 7B), implying that the five ZNF TFs indeed involve in the anti-SARS-COV-2 infection. Importantly, we found that target genes of females could be enriched in some specific items, such as ‘myeloid cell differentiation’, ‘myeloid lymphocyte differentiation’, ‘interaction with host’ and ‘negative regulation of virus release’ (Figure 7A). In contrast, these target genes of males are mainly enriched in the ‘positive regulation cytokine production’, ‘cellular response to lipopolysaccharide’ and ‘cellular response to biotic stimulus’ (Figure 7B). Especially, some target genes of females can be enriched in some certain molecular function, such as helicase activities, ‘histone binding’ and ‘ubiquitin-like protein ligase binding’ (Figure 7B).

BALF scRNA-seq reveals the anti-SARS-CoV-2 role of ZNF proteins

The 60 326 BALF single cells could be further clustered into 27 clusters (Figure S10, see Supplementary Data available online at http://bib.oxfordjournals.org/). According to the specific marker genes, 27 cell clusters were identified as 12 main cell types (Figures 8A and B and S11, see Supplementary Data available online at http://bib.oxfordjournals.org/). We found that the overall ZNF protein activity of different cell types undergone extensive changes after SARS-CoV-2 infection (Figure 8C). Compared with healthy people, the ZNF protein activity of the lung epithelial cells was significantly decreased in mild and severe patients (Figure 8C). Conversely, the ZNF activity of some immune cells, such as macrophages, CD8 T cells, CD4 T cells and NK cells, is, respectively, significantly increased in mild patients, but significantly reduced in severe patients (Figure 8C). Especially, the expression level and proportion of ZNFX1 in mild and severe patients were significantly increased, and similar expression trends also existed in other key immune genes (Figure 8D). The overall ZNF protein activity in the lung immune cells of mild and severe females was also higher than males, especially macrophages, T cells and NK cells (Figure 8E and F).

Analysis of clinical indicators of male and female patients

Herein, we analyzed the TBNK cells data (T cells, B cells and NK cells) from 548 COVID-19 patients. Fascinatingly, female patients have higher immune abundance of CD 3+, CD4+ T cells, B cells and CD4/CD8 ratio than males (Figure 9A and E), but males have more NK cells than females (Figure 9F). Female patients have more naïve B cells and effectors CD8 T cells than males (Figure 9G and H). Although no significant differences existed in other 25 cell subsets (Figures S12–S14, see Supplementary Data available online at http://bib.oxfordjournals.org/), female patients do have higher proportion and stronger activity of lymphocytes than males.

Herein, we further found that the incidence of critical illness and death patients is significantly increased with abnormal decrease of lymphocyte cell abundances (hazard ratio: 2.7330), indicating that lymphocytes could serve as a good protective factor for COVID-19 patients (Figure S15, see Supplementary Data available online at http://bib.oxfordjournals.org/). The lymphocyte abundance of male patients is always lower than females in earlier multiple tests, but the lymphocyte abundance gradually tends to be consistent between them with the treatment progresses (Figure 9I). Remarkably, lymphopenia is more detrimental to male COVID-19 patients than females with hazard ratio (male: 3.2977 versus female: 1.9776) (Figure S15, see Supplementary Data available online at http://bib.oxfordjournals.org/). Finally, we analyzed the effect of major lymphocyte subtypes on the prognosis of COVID-19 patients. Our results revealed that CD3, CD4, T cells and B cells as protective factors were significantly associated with the prognosis of COVID-19 patients, but CD8 T cells, NK cells and CD4/CD8 ratio were no significant effect on the prognosis of COVID-19 patients (Figure 9J).