Pyroptotic Patterns in Blood Leukocytes Predict Disease Severity and Outcome in COVID-19 Patients

The global coronavirus disease 2019 (COVID-19) pandemic has lasted for over 2 years now and has already caused millions of deaths. In COVID-19, leukocyte pyroptosis has been previously associated with both beneficial and detrimental effects, so its role in the development of this disease remains controversial. Using transcriptomic data (GSE157103) of blood leukocytes from 126 acute respiratory distress syndrome patients (ARDS) with or without COVID-19, we found that COVID-19 patients present with enhanced leukocyte pyroptosis. Based on unsupervised clustering, we divided 100 COVID-19 patients into two clusters (PYRcluster1 and PYRcluster2) according to the expression of 35 pyroptosis-related genes. The results revealed distinct pyroptotic patterns associated with different leukocytes in these PYRclusters. PYRcluster1 patients were in a hyperinflammatory state and had a worse prognosis than PYRcluster2 patients. The hyperinflammation of PYRcluster1 was validated by the results of gene set enrichment analysis (GSEA) of proteomic data (MSV000085703). These differences in pyroptosis between the two PYRclusters were confirmed by the PYRscore. To improve the clinical treatment of COVID-19 patients, we used least absolute shrinkage and selection operator (LASSO) regression to construct a prognostic model based on differentially expressed genes between PYRclusters (PYRsafescore), which can be applied as an effective prognosis tool. Lastly, we explored the upstream transcription factors of different pyroptotic patterns, thereby identifying 112 compounds with potential therapeutic value in public databases.

A possible mechanism linking cytokine storm to organ damage is inflammatory cell death, namely pyroptosis and necroptosis. Pyroptosis has been intensely studied recently. Some patients with severe COVID-19 may develop a systemic cytokine storm because SARS-CoV-2 promotes cytokine storms by inducing pyroptosis in pro-inflammatory blood-born immune cells (16)(17)(18)(19). However, only a few studies about necroptosis and cytokine storm in COVID-19 have been published thus far.
Pyroptosis is a mechanism of programmed cell death characterized by the inflow of sodium ions and water mediated by gasdermin proteins, resulting in cell membrane rupture, excessive cell swelling, and spontaneous release of cytosolic contents into the extracellular space (20). Gasdermin proteins, which consist of an N-terminus with membrane pore-forming activity and an inhibitory C-terminus, are the key regulators of pyroptosis. Upon inflammasome activation, caspase proteins, including caspase-1 and other non-canonical inflammasome caspases (e.g., caspase-4, caspase-5, and caspase-11), cleave gasdermin into two parts (21), thereby unleashing the poreforming activity of the N-terminus. This N-terminus fragment of gasdermin binds to the cell membrane, forming pores and leading to pyroptosis (22). Pyroptosis triggers the rapid release of a slew of alarmins including, cytokines (IL-1b, IL-18), chemokines, and damage-associated molecular patterns (DAMPs), prompting an immediate response from surrounding immune cells and triggering a pyroptotic chain reaction (23). Thus, pyroptosis plays a key role in the emergence of a cytokine storm, according to recent research (17,19,24).
Although pyroptosis is crucial for innate immunity (25,26), extensive pyroptosis can cause tissue inflammation, organ failure and death, as found in various diseases (27). For example, in atherosclerosis, cholesterol crystals and oxidized low-density lipoprotein cause macrophage pyroptosis, which leads to a massive release of cytokines, promoting inflammation and disease progression (28). In line with these findings, NLRP3 inflammasome or ASC inhibition, which prevent macrophage pyroptosis, can lower infarct size and improve heart function in an animal model of myocardial infarction (29,30). Similar results have also been observed in alcoholic hepatitis (31), lupus erythematosus (32,33), and even in the central nervous system (34).
In COVID-19 patients, various cells undergo pyroptosis, including leukocytes (monocytes, macrophages, mucosal-associated invariant T cells) and other type of cells (adipocytes, lung epithelial cells and endothelial cells) (16,(35)(36)(37)(38)(39). On the one hand, SARS-CoV-2 nucleocapsid can prevent Gasdermin D cleavage, thus reducing host pyroptosis and suppressing the immune response (40), in addition to inhibiting coronavirus infection by promoting non-classical secretion of b-interferon (41). On the other hand, SARS-COV-2 can stimulate macrophage GSDMD-mediated pyroptosis, which leads to the rapid release of pro-inflammatory cytokines and to a cytokine storm (17). In addition, the synergistic effect of TNF-a and IFNg can trigger GSDMD-mediated pyroptosis and promote a cytokine storm, thereby increasing mortality among COVID19 patients (42). However, in cats and dogs, deficiencies of the inflammasome and pyroptosis pathways (cats and tigers do not express AIM2 and NLRP1, and dogs do not express AIM2 and have a shorter form of NLRC4 than humans) may provide an evolutionary advantage against SARS-CoV-2 by reducing cytokine storm-induced host damage (43). Therefore, the role of pyroptosis in COVID-19 remains complex, requiring more comprehensive studies.
Based on transcriptome data of patients with or without COVID-19 available in public databases (GSE157103), we found that the leukocytes of ARDS patients with COVID-19 have considerably higher pyroptotic markers than patients without COVID-19. Moreover, at least two different patterns of pyroptosis occur in patients with COVID-19, one correlated with a poor prognosis and the other with a benign prognosis. These two pyroptosis patterns may be regulated by different upstream transcription factor networks, which could prove therapeutically valuable for drug development.

MATERIALS AND METHODS
Obtaining RNA-seq Data from the GEO Dataset From the GEO dataset, we retrieved RNA-seq data (GSE157103) of 126 ARDS patients, namely 100 COVID-19 patients and 26 non-COVID-19 patients, in addition to their clinical data, including gender, age, underlying disease status (diabetes), coagulation (D-dimer, ferritin, CRP, procalcitonin, fibrinogen), and hospital-free days post 45-day follow-up (HFD45), among other parameters. More specifically, HFD45 is defined as the number of days patients lived outside of a hospital from enrollment through death or the end of follow-up (44). The higher HDF45 is, the milder the disease and the better the prognosis will be. proteins were sorted from large to small by log2FC (not absolute value). Then, using the "org.Hs.eg.db" and "clusterProfiler" packages, gene set enrichment analysis (GSEA) was performed based on MSigDB gene sets C2, C5 and C7. Significant gene sets were identified when |Normalized Enrichment Score (NES)|>1 and False Discovery Rate (FDR) < 0.25.

Function and Pathway Analysis of DE Immune Genes
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the "org.Hs.eg.db" and "clusterProfiler" packages. GO terms and KEGG terms were identified as significantly enriched when p.adjust < 0.05.

Estimation of Immune Cell Infiltration Fractions
Single-sample gene set enrichment analysis (ssGSEA) and "Cibersort" were used to analyze the immune cell infiltration fractions. The former was based on the list of Pan-cancer Immune Metagenes (45,46).

Construction of the PYRscore
The differentially expressed genes (DEGs) of two PYRclusters were identified using the "limma" package. Using the median of HFD45 as the cutoff value, COVID-19 patients were divided into two groups. The DEGs in these two groups were calculated using the limma package and subsequently applied for PCA analysis. PC1 and PC2 were used to construct the PYRscore (47).

Construction of the PYRsafescore Model
Based on the log2(TPM) of the 570 DEGs between PYRclusters and the HDF45 of each patient, we used the "glmnet" package to build a PYRsafescore model by LASSO regression. We determine the signatures of the model by selecting the lambda value with the smallest mean-squared error by 20-fold-cross-validation. The coefficients of the final signatures were used to calculate the PYRsafescore as follows: protective score = ∑ Coefficient i × Expression level of signature i . Using the "caret" package, 100 patients were randomly divided into a training group and a test group with a ratio of 2:1. The model built with the training group data was validated in the test group. We used "ROCR" packages to plot receiver operating characteristic (ROC) curves and to calculate area under the curve (AUC) scores to evaluate model performance.

Transcription Factors Enrichment
First, based on the transcription factor targets (TFT) and their gene sets in MSigDB (48), we used the "clusterProfiler" package to enrich the transcription factors from DEGs between PYRclusters. We then used the "CoRegNet" package to enrich the transcription factor co-regulatory network in COVID-19 patients from PYRcluster1 and PYRcluster2 based on the dataset of transcription factor targets (CHEA, ENCODE, and JASPAR Predicted Transcription Factor Targets, MotifMap and TRANSFAC Predicted Transcription Factor Targets, and TRANSFAC Curated Transcription Factor Targets) from Harmonizome (49,50). Lastly, the network graph was plotted using the "ggraph" package.

Search for Drugs Targeting Transcription Factors
We used the transcription factors that we screened as keywords to search for the corresponding compounds on ChEMBL (51), thus screening active or repressed compounds based on their role in pyroptosis.

Statistical analysis
The Wilcoxon sum-rank test and the t-test were used to compare different groups, and the Pearson's product-moment correlation test was used for correlation analysis. All statistical tests were two-sided, and a significant difference was defined as a p-value of 0.05. Power calculations were performed using the following R packages: "pwr" and "rstatix" at sig.level=0.05.

Transcriptome Data Reveal the Pyroptosis Characteristic of Blood Leukocytes from COVID-19 Patients
A GEO dataset provided RNA-seq data and clinical data from 126 samples of 100 patients with COVID-19 and 26 patients without COVID-19 (GSE157103) (44). Initially, we assessed the expression levels of pyroptosis-related gene sets of all patients based on prior studies (52).
By gene set variation analysis (GSVA), we studied changes in the biological function of leukocytes between the two types of patients. Mismatch repair, homologous recombination, replication, cell cycle, and p53 signaling are more enriched in leukocytes of COVID-19 patients, indicating severe cell damage during viral infection and ongoing damage repair ( Figure 2A and Supplementary Table 1).
By single-sample gene set enrichment analysis (ssGSEA), we also compared the proportions of 28 immune cell types between the two groups and found that the numbers of various immune cells are much higher in COVID-19 patients than in non-COVID-19 patients, indicating a highly active immune response in COVID-19 patients (46) ( Figure 2B and Supplementary Table 2). Surprisingly, the numbers of some immune cells, such as macrophage and cd56 dim natural killer cells, were lower in COVID-19 patients than in the control group. This difference could be due to cell death caused by massive viral infection and increased pyroptosis (17,55). Human dendritic cells and T lymphocytes (including CD4 + and CD8 + ) undergo pyroptosis via the AIM2-Caspase1-gasdermin D and the CARD8-Caspase1-gasdermin D axes (56,57), respectively. Pyroptosis has also been identified in macrophages and neutrophils (58,59). Although NK cells have not been documented to undergo pyroptosis on their own, they can participate in this process, playing a key role (60). These results suggests that blood leukocytes either directly undergo pyroptosis or have a synergistic role with pyroptosis in COVID-19 patients.

COVID-19 Patients Showed Different Patterns of Pyroptosis
By non-negative matrix factorization based on 35 pyroptosisrelated genes, we clustered the 100 COVID-19 patients into two clusters, PYRcluster1 and PYRcluster2 ( Figure 3A and Supplementary Figure 1A). The best clustering result was found when k = 2, with no differences in age, gender, days admitted before enrollment, replacement therapy (preenrollment) or underlying disorders between the two clusters, which we called PYRclusters (Supplementary Figures 1B-E, 3A, B).
Subsequently, we found that PYRcluster2 has a higher hospital-free days post 45-day follow-up (HFD45), ventilatorfree days and a lower proportion of mechanical ventilation than PYRcluster1 ( Figures 3D, E and Supplementary Figure 3C  which suggests a better prognosis. The blood D-dimer level of PYRcluster1 is significantly higher than that of PYRcluster2 ( Figure 3F), indicating hypercoagulation. Although albumin and hemoglobin of the two PYRclusters are mostly below the normal range (green dashed line), PYRcluster1 deviates further from the normal range than PYRcluster2. Other clinical features do not differ significantly between the two clusters ( Supplementary Figures 2A-C). Based on these results, the highly expressed pyroptosis-related genes of PYRcluster1 may be associated with a poor prognosis. In line with our results, AIM2 and NLRC4 deficiency in dogs and cats provide a protective effect against SARS-CoV-2 by reducing cytokine storm-induced host damage (43). We further determined differentially expressed genes (DEGs) using the "limma" package, with cutoff criteria of | logFC| >1 and p=0.05, totaling 570 DEGs, and employed Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichments ( Figure 4A and Supplementary Figure 2E). In addition, we determined the logFC of 736 proteins between two PYRclusters (FC=PYRcluster2/ PYRcluster1) using "limma" from the proteomic data from the same study as the transcriptome data. Subsequently, we employed gene set enrichment analysis (GSEA) based on the log2FC of these proteins (Supplementary Figure 2D, Supplementary Table 4).
Immune responses such as antigen recognition/presentation, immune cell activation, migration, and replication are relatively enhanced in PYRcluster2 ( Figure 4A and Supplementary Figure 2E). Both transcriptome and proteomic results revealed hyperinflammation in PYRcluster1 (Supplementary Figure 2D, Supplementary Table 6). KEGG pathway enrichment data showed that "Coronavirus disease-COVID-19" pathways are markedly upregulated in PYRcluster2, underscoring a highly activated immune process unique to PYRcluster2. PYRcluster1 has a higher NO synthesis level than PYRcluster2, which may be related to the antiviral capacity of its patients (61)(62)(63). The highly activated stress response pathway may associated with severe damage caused by the virus and lead to higher levels of immune cell apoptosis in PYRcluster1 (Supplementary Figure 2D, Supplementary Table 6). Furthermore, PYRcluster1 has a markedly increased expression of cytokines, including IL-1, IL-6, IL-8, IL-10, and TNF (Supplementary Table 6), which are typical components of a cytokine storm and result in a poor prognosis (64,65). Consistent with an increased D-dimer ( Figure 3F), the neutrophil extracellular traps (NETs) that promote blood coagulation are highly expressed in PYRcluster1, promoting venous thrombosis and leading to poor prognosis (66,67). In addition, AIM2, CXCR2 and UBE2W, which were significantly highly expressed in PYRcluster1, were included in the 20 genes associated with distinctly methylated CpG sites between mild and severe COVID-19 patients (68). Their odds ratios were all greater than 1, indicating that their downregulation is beneficial to COVID-19 patients. In conclusion, the two PYRclusters of COVID-19 patients exhibited distinct pyroptotic patterns and clinical features.

Different Pyroptotic Patterns are Associated with Different Leukocytes and Opposite Prognosis
SARS-CoV-2 viruses infect leukocytes and lead to immunodeficiency (69)(70)(71). We speculate that leukocyte pyroptosis helps to destroy the virus protective niche and release viruses from cells, thereby enhancing viral clearance and immune recovery. As a result, viruses released by pyroptosis may be further removed by phagocytic cells via phagocytosis.
We first assessed the proportion of different immune cells in the two PYRclusters by ssGSEA. Immune cells highly associated with antivirals, such as activated B, CD4/8 + T, Treg, and NK cells are highly expressed in PYRcluster2 ( Figure 4B and Supplementary Figure 4A, Supplementary Table 7). PYRcluster1 had a higher proportion of pro-inflammatory neutrophils. Using "Cibersort", we discovered PYRcluster2 had a higher proportion of antiinflammatory M2 macrophages (Supplementary Figure 4A). The hemogram percentages of several cells, such as neutrophils, lymphocytes and monocytes, were consistent with the results of ssGSEA (Supplementary Figure 3F). We then used the "estimate" package to score the two clusters and found that PYRcluster2 shows a greater increase in immune cell infiltration than PYRcluster1 ( Figure 4C). As shown in Figure 4D, the expression of characteristic pyroptosis-related genes of PYRcluster1 were significantly positively correlated with its expression of characteristic immune cells. For example, PYCARD, GPX4, GSDMD, GZMA, TNF, NOD1, IL6, NLRP7, CASP6, GSDME, PJVK, GSDMA, and NLPR2 are positively correlated with NK, NKT, regulatory T, gd T, activated B, and Th17 cells, MDSCs, and monocytes ( Figure 4D). Both these pyroptosis-related genes and immune cells are highly expressed in PYRcluster2 (Figures 3C, 4B). Furthermore, leukocytes have stronger phagocytic activity in GO enrichment results in PYRcluster2 than in PYRcluster1 (Supplementary Table 6). In contrast, pyroptosis in PYRcluster1 may produce several pathogenand damage-related molecular patterns that increase cytokine storm, leading to multiple organ failure and poor prognosis.
We subsequently performed unsupervised clustering of all COVID-19 patients using 570 DEGs and all genes, respectively, yielding two more clusters: DEG and All-Gene clusters. The heatmap demonstrates that these additional clusters match a previous clustering based on the 35 pyroptosis-related genes ( Figure 5A and Supplementary Figure 4B-D). These data suggest that distinct patterns of pyroptosis occur in COVID-19 patients, which can be represented by 35 pyroptosis-related genes.
To better elucidate differences in pyroptotic patterns between PYRcluster1 and PYRcluster2 and their correlation with prognosis, we created a pyroptosis score (pyrscore) (Supplementary Figure 4E). As shown in Figure 5B, PYRcluster1 has a higher score than PYRcluster2. HFD45 and ventilator-free days are negatively correlated with pyroptosis scores (Figures 5C, D), whereas sofa, APACHE-II, D-dimer, and CRP levels are positively correlated ( Supplementary  Figures 5A-D). In conclusion, PYRcluster1 and PYRcluster2 have different levels of immune response to SARS-COV-2 and pyroptotic patterns, resulting in distinct prognoses.

Development of A Predictive Pyroptotic Prognosis Model
Given that different pyroptotic patterns may have a significant impact on the prognosis of COVID-19 patients, we created a PYRsafescore model based on the HFD45 of COVID-19 patients and DEGs across PYRclusters by least absolute shrinkage and selection operator (LASSO) regression analysis. In a ratio of 2:1, 100 COVID-19 patients were divided into a training group and a test group, and the model was obtained in the training group. "Lambda-min" was chosen as the best value in the cross-validation procedure ( Figure 5E and Supplementary Figures 5E). Lastly, based on the log2 value of the expression level of 10 genes, we established the following scoring model: The area under the ROC curve (AUC) values of the model in the training and test groups were 0.907 and 0.879, respectively, indicating that our predictive model performs well ( Figure 5F). HFD45 and PYRsafescore are positively correlated ( Figure 6A), which suggests that a higher PYRsafescore indicates a better prognosis. In addition to HDF45, the correlation of other clinical variables with PYRsafescore, including sofa, ventilator-free days, APACHE-II, and CRP levels, also demonstrates that our model works effectively (Figures 6B, C and Supplementary  Figures 6A, B). Additionally, the expression of these 10 genes between the two PYRclusters also has significant differences ( Figure 6D). KEGG analysis indicates that they were closely related to the patient's immune response ( Figure 6E). Frontiers in Immunology | www.frontiersin.org July 2022 | Volume 13 | Article 888661 In conclusion, our results proved that the newly created prognosis model has a considerable clinical predictive value.

Transcriptional Regulatory Networks and Potential Drugs in Different Pyroptosis Patterns
Using the "clusterProfiler" package, we first enriched DEGs across PYRclusters for transcription factors based on MSigDB Collections: "regulatory target gene sets" (Figure 6F). We further enriched the transcription factor regulatory networks for all transcriptome data of the two PYRclusters using the "CoRegNet" package. PYRcluster1 has the regulatory network with MNDA, TSC22D3, HMGB2, FOS, EEF1A1, TRIM22, NFKBIA as transcription factors, and PYRcluster2 has the regulatory network with FOS, MNDA, EEF1A1, DAZAP2, DDIT3, HCLS1, NFKBIA, TSC22D3, PTMA, and TRIM22 ( Figure 7A and Supplementary Figure 6C  Once SARS-COV-2 infects cells, EEF1A1 is critical for viral replication. Drug targeting EEF1A has robust antiviral effects in vitro (72). The nuclear factor-kappaB (NF-kappaB)/REL family of transcription factors are activated in response to DNA damage to regulate inflammation and apoptosis resistance (73,74). Since NFKBIA is highly expressed in COVID-19 patients (75) and tripartite motif-containing (TRIM) 22 can activate NF-kappaB (76) to protect the host from viral infection (77), these two genes have a huge impact on immune disparities between the two PYRclusters. MNDA is a member of the family of hematopoietic interferon (IFN)-inducible nuclear proteins that promotes the degradation of the anti-apoptotic factor MCL-1 and apoptosis in myeloid cells (78,79).
We then observed the correlation between these transcription factors and the expression of differentially expressed pyroptosis-related genes between PYRclusters and the clinical characteristics of all patients. The correlations between transcription factors and pyroptosis-related genes vary greatly, depending on the PYRclusters ( Figure 7B Frontiers in Immunology | www.frontiersin.org July 2022 | Volume 13 | Article 888661 related genes in PYRcluster2 also have a positive effect on prognosis ( Figure 7C). Taken together, these findings show that distinct pyroptotic patterns may result from different upstream transcriptional regulation pathways. In light of this, we searched the ChEMBL database for compounds that promote or inhibit appropriate transcription factors based on diverse roles in prognosis (51). We screened a total of 112 compounds (Supplementary Table 9) and found that CHEMBL348436 (also known as Cirsimaritin) has the potential to regulate blood leukocyte pyroptosis in COVID-19 patients while simultaneously improving prognosis through FOS and NFKB1A inhibition. In fact, drugs targeting FOS have therapeutic effects in COVID-19 patients (80).

DISCUSSION
Pyroptosis, a mechanism of programed cell death which leads to cell swelling and lysis, plays a key role in innate immunity by disrupting the pathogen replication niche and killing intracellular bacteria through pore-induced intracellular traps (25,26). However, excessive pyroptosis may trigger an overactive inflammatory response, resulting in a cytokine storm and severe organ damage through IL-6, TNF and NETs (81)(82)(83).
The occurrence of a cytokine storm is a major factor in the progression of moderate-to-severe COVID-19. In the therapy of COVID-19, multi-organ failure induced by cytokine storm has become a significant issue (84). Despite the relevance of pyroptosis in the treatment of severe COVID-19 patients, research on COVID-19 and pyroptosis is currently limited. Some studies suggest that the elevated pyroptosis is not conducive to the treatment of the disease but closely related to SARS-CoV-2 infection and cytokine storm (17,35,42,43). By contrast, other studies show that pyroptosis can also be beneficial in fighting SARS-CoV-2 infection (40,41). A dual role for NLRP3 was reported in a recent study according to which inflammasome-dependent pyroptosis contributes to the hyperinflammatory state of the lungs. However, pyroptosis can release infectious virus, preventing a productive viral cycle, which can help to eliminate viruses (85). In conclusion, the role of pyroptosis in COVID-19 remains unclear.
In combination with other clinical information, we found that PYRcluster2 patients had a better prognosis, including a longer HFD45 and a lower ICU hospitalization rate. PYRcluster2 patients also had more immune cells and a higher immune score. Furthermore, the expression of immune cells was highly correlated with the expression of pyroptosis-related genes. To better elucidate the pyroptotic patterns and prognosis, we calculated the "pyrscore" to characterize different pyroptotic patterns and the "PYRsafescore" to better predict prognosis and assist clinical treatment. Higher PYRsafescore scores mean a better prognosis. Lastly, by transcription factor enrichment, we identified the upstream transcription factors that regulate different pyroptotic patterns and screened a series of compounds with therapeutic potential in public databases.
Currently, only a few drugs are available for controlling SARS-CoV-2 infection, including monoclonal antibodies that neutralize viral proteins (90)(91)(92), drugs that inhibit viral replication (93), and new oral drugs from Pfizer and Merck, namely PAXLOVID and Molnupiravir, which inhibit viral replication and viral proteases, respectively (94,95). Methods for managing excessive inflammatory response, which is the leading cause of severe COVID-19, are very limited and less effective (96): For example, heparin is widely used to prevent blood clots, and some immunosuppressants such as dexamethasone, IL-6 monoclonal antibodies and JAK kinase family inhibitors are used to inhibit inflammation (12)(13)(14)(15)(97)(98)(99). Our results suggest that pyroptosis plays a key role in the generation of a hyperinflammatory immune response. Therefore, therapeutic strategies targeting pyroptosis have potential value in managing inflammation and hence reducing COVID-19 severity and mortality.
In conclusion, our data reveal different patterns of pyroptosis of blood leukocytes in patients with COVID-19, which is closely related to their prognosis. We speculate the mechanism underlying diverse prognoses is shown in Figure 8. Prognosis prediction models developed based on different pyroptosis patterns are highly valuable for COVID-19 treatment. In addition, compounds that target the transcription factor network that regulates the pyroptotic process may help to develop new drugs for the treatment of patients with severe COVID-19.

AUTHOR CONTRIBUTIONS
YT, PZ, and LC designed and supervised the study, PZ collected the data, PZ performed the first part of the analysis, and YT completed the analysis and drafted the manuscript. YT, PZ, and QL prepared the original draft and wrote, reviewed, and edited the manuscript. JX and LC revised the manuscript and provided analytical technical support. All authors contributed to the article and approved the submitted version.