Predicting the Prognosis of Esophageal Adenocarcinoma by a Pyroptosis-Related Gene Signature

Esophageal adenocarcinoma (EAC) is a highly malignant type of digestive tract cancers with a poor prognosis despite therapeutic advances. Pyroptosis is an inflammatory form of programmed cell death, whereas the role of pyroptosis in EAC remains largely unknown. Herein, we identified a pyroptosis-related five-gene signature that was significantly correlated with the survival of EAC patients in The Cancer Genome Atlas (TCGA) cohort and an independent validation dataset. In addition, a nomogram based on the signature was constructed with novel prognostic values. Moreover, the downregulation of GSDMB within the signature is notably correlated with enhanced DNA methylation. The pyroptosis-related signature might be related to the immune response and regulation of the tumor microenvironment. Several inhibitors including GDC-0879 and PD-0325901 are promising in reversing the altered differentially expressed genes in high-risk patients. Our findings provide insights into the involvement of pyroptosis in EAC progression and are promising in the risk assessment as well as the prognosis for EAC patients in clinical practice.


INTRODUCTION
Esophageal cancer is one of the most common malignancies worldwide, accounting for approximately 604,100 new cases and 544,076 deaths per year over the world (Sung et al., 2021). Esophageal adenocarcinoma (EAC) and esophageal squamous cell carcinoma (ESCC) composite the principle histologic subtypes of esophageal cancer, in which the incidence of EAC in western countries has increased dramatically in the last decades (Klingelhöfer et al., 2019). Despite therapeutic advances in surgery, radiotherapy, chemotherapy, and targeted drugs, the 5-year survival of esophageal cancer remains less than 20% (Alsop and Sharma, 2016). In consequence, biomarkers and effective models are urgently needed to predict the prognosis of EAC and provide insights into targeted therapy.
Pyroptosis is a proinflammatory form of regulated cell death, relying on the enzymatic activity of inflammatory proteases that belong to the caspase family (Vande Walle and Lamkanfi, 2016). Pyroptosis is featured with swift plasma-membrane rupture and subsequent release of proinflammatory intracellular contents, which is distinct from apoptosis (Bergsbaken et al., 2009). Studies evaluating the role of pyroptosis in neurological, infectious, autoimmune, cardiovascular, and oncologic disorders have been emerging in recent years (Yu et al., 2021). Activation of the canonical inflammasome pathway is the basis of pyroptosis, in which patternrecognition receptors (PRRs), for example, Toll-like receptors (TLRs), nucleotide-binding oligomerization domain-like receptors (NLRs), and absent in melanoma 2 like-receptors (ALRs) recognize pathogen-associated molecular patterns (PAMPs) or nonpathogenrelated damage-associated molecular patterns (DAMPs) to activate inflammasomes and facilitate caspase-1 activation (Xia et al., 2019). Direct activation of caspase-4/5/11 under lipopolysaccharide (LPS) is involved in the noncanonical pyroptosis pathway, which is independent of the inflammasome complex (Shi et al., 2014). The gasdermin (GSDM) family proteins serve as the main mediators of pyroptosis, which are proteolytically activated by proteases and induce the formation of plasma membrane pores, leading to cell swelling and lysis (Van Opdenbosch and Lamkanfi, 2019;Tsuchiya, 2020). Due to the pivotal role of GSDM family proteins, pyroptosis is defined by some researchers as gasdermin-mediated programmed cell death (Shi et al., 2017).
However, despite the fact that research is emerging in ESCC, the role of pyroptosis in esophageal cancer remains largely unknown, and none of the previous publications have comprehensively evaluated the pyroptosis-related genes in EAC. Therefore, we performed a comprehensive evaluation of pyroptosis-related genes in EAC, in order to develop a pyroptosis-gene-based modality to predict the prognosis of the patients, and provide insights into the correlations between pyroptosis and tumor immune microenvironment.

Datasets
The RNA-sequencing (RNA-seq) data of 87 patients (78 with EAC; 9 normal samples) and the corresponding clinical information from The Cancer Genome Atlas (TCGA) database were retrieved on May 20, 2021 (https://portal.gdc.cancer.gov/ repository). The DNA microarray and clinical features of the validation cohort were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/, ID: GSE13898). The initial inclusion criteria were as follows: 1) patients with EAC; 2) patients with clear data for overall survival and survival status; and 3) patients with available gene expression data. The exclusion criteria were as follows: 1) patients with ESCC; 2) patients with incomplete data for overall survival, survival status; and 3) patients without gene expression data. As described in the following sections, further analysis based on clinicopathological characteristics was performed in patients with complete clinical data including age, gender, and stage. Patients with survival time of less than 30 days were excluded.
Yin, 2017; Karki and Kanneganti, 2019;Xia et al., 2019). The microarray data from the GSE13898 cohort were normalized using the quantile normalization method, and the expression levels of genes were transformed to a log 2 base for further analysis (Kim et al., 2010). The package "limma" was used to explore DEGs with the threshold of p value <0.05 (Ritchie et al., 2015). Probes with missing information for gene expressions in >20% samples were removed. The correlations of selected genes were evaluated by the "ggcorrplot" package (Kassambara and Kassambara, 2019). Protein-protein interaction (PPI) networks were created by Search Tool for the Retrieval of Interacting Genes (STRING) and the "igraph" package (Csardi and Nepusz, 2006;Szklarczyk et al., 2019).

Development and Validation of the Pyroptosis-Correlated Gene Prediction Model for Prognosis
Cox regression analysis was employed to evaluate the value of pyroptosis-related genes for prognosis. The DEGs were identified for further analysis. The LASSO Cox regression analysis was employed to construct a refined model for prognosis using the R package "glmnet" (Friedman et al., 2010). The calculation of the risk score was performed using the following formula: risk score i (n 1) Coef ipXi (Coef i indicates the coefficient, and Xi indicates the gene expression levels after standardization). The EAC patients were classified into low-and high-risk groups based on the median risk score, and Kaplan-Meier analysis was used to compare the overall survival (OS) between the two groups. Principal component analysis (PCA) was used to assess the separability of the two groups by the "prcomp" function. The R packages "survival," "survminer," "timeROC," and "riskRegression" were utilized for receiver operating characteristic (ROC) curve graphing and area under curve (AUC) calculation for 1, 2, and 5 years (Blanche et al., 2013;Therneau and Lumley, 2015;Kassambara et al., 2017;Ozenne et al., 2017). A nomogram model with clinical features including stage and risk score was constructed by the R packages "rms," "foreign," and "survival" (Therneau and Lumley, 2015;Harrell et al., 2017;Team et al., 2020). The calibration curve and detrended correspondence analysis (DCA) were performed using the "rms" package (Harrell et al., 2017). An EAC cohort (GSE13898) from the GEO database was used for validation, and the risk score was calculated by the same methods described above to divide the cohort into two subgroups (low risk and high risk).

Prognostic Analysis of the Variables
Clinical data (age, gender, and stage) were extracted from patients in the TCGA and GSE13898 cohorts. The clinicopathological characteristics of EAC patients with complete data for further analysis were described in Supplementary Table S2. Variables including gender, stage, and risk score were analyzed in the regression model by univariate and multivariate Cox regression analysis.

Methylation Analysis
For the genes included in the signature, the cBio Cancer Genomics Portal (cBioPortal) database (http://www.cbioportal. org/) was used for exploring the correlation between methylation alterations and gene expressions in the TCGA esophageal adenocarcinoma cohort. The MEXPRESS database (http:// mexpress.be/) was utilized for further assessment of the correlation between the precise genomic location of DNA methylation and altered levels of gene expression. p < 0.05 and R 2 > 0.25 were considered as significant correlation.

Tumor Microenvironment Analysis
The Tumor Immune Estimation Resource (TIMER) database (https://cistrome.shinyapps.io/timer/) was utilized to assess the correlation between tumor-infiltrating immune cells and expressions of selected genes (Li et al., 2017). Estimation Resource (TIMER) was used to compare the immune scores of the four subtypes. The CIBERSORT algorithm was used to further explore the composition and differences in the fraction of 22 immune cell types between two subgroups classified by risk scores (Chen et al., 2018).

Enrichment Analysis
Patients with EAC in the TCGA cohort were divided into two groups based on the median risk score. The DEGs between the low-and high-risk groups were extracted by |log2FC| ≥ 1 and p value <0.05. GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was performed by the R package Frontiers in Pharmacology | www.frontiersin.org November 2021 | Volume 12 | Article 767187 "clusterProfiler," and the results were visualized using the "GOplot" package (Yu et al., 2012;Walter et al., 2015).

Connectivity Map Analysis
CMap analysis was performed and visualized in https://clue.io/ (Lamb et al., 2006), The top 150 upregulated and downregulated genes were selected according to the |logFC| values of DEGs for CMap analysis to identify a shortlist of drugs. According to the pattern-matching algorithms, positive scores indicate the induction effect of the small molecules on the signature, while negative scores indicate the inhibition effect. The drugs were further selected based on the negative scores.

Statistical Analysis
Statistical analyses were performed by R (version 4.1.0). Student's t-test was applied to compare the differences in gene expression between tumor and normal tissues, while categorical variables were compared using Pearson's chi-square test. The OS of patients between low-and high-risk groups were compared by the Kaplan-Meier method with log-rank test. The Cox regression analysis was performed to evaluate the independent prognostic factors for survival. The Wilcoxon test was used to compare the immune cell infiltration between groups.

Code Availability
The R code used in this study is available from the corresponding author upon reasonable request.

Identification of DEGs Between EAC and Normal Tissues
The expression levels of 58 pyroptosis-correlated genes were examined in the TCGA data of 78 EAC and 9 normal tissues. Ten DEGs were identified (|log2FC| ≥ 1 and p value <0.05), and all of them (CASP1, CASP5, GSDMB, GZMB, IL1B, NLRP6, PYCARD, TNF, TREM2, and ZBP1) were upregulated in the tumor group. The expression profiles of DEGs were demonstrated in Figure 1A (red color represents a higher expression level; blue color represents a lower expression level). Figure 1B showed the correlation network of DEGs in the TCGA data, indicating that GSDMB expressions are strongly correlated with CASP1 (r 0.80, p < 0.05) and TREM2 (r 0.67, p < 0.05) expressions. In addition, the expression of CASP1 is significantly correlated with that of CASP5 (r 0.64, p < 0.05). The PPIs of DEGs were presented in Figure 1C, in which the interaction score was set as 0.4. The correlation between CASP1 and CASP5 was consistent in the protein level.

Construction of Prognostic Model Based on DEGs
A total of 65 EAC patients with available survival data were included in our study. Univariate Cox regression analysis was initially performed to assess the prognostic value of DEGs Frontiers in Pharmacology | www.frontiersin.org November 2021 | Volume 12 | Article 767187 5 ( Figure 2A). Among them, six genes (CASP1, CASP5, GSDMB, IL1B, PYCARD, and ZBP1) were with p value <0.2, and higher expressions of CASP1, CASP5, and IL1B were associated with increased risk (HR > 1), while upregulated expressions of GSDMB, PYCARD, and ZBP1 were correlated with lower risk (HR < 1). Subsequently, LASSO Cox regression analysis retrieved five genes for prognostic model construction based on the optimum λ value ( Figures 2B,C). The calculation of the risk score was as follows: Risk score (0.042 × expCASP1)＋(−0.025 × expGSDMB) + (0.021 × expIL1B) + (−0.037 × expPYCARD) + (−0.243 × expZBP1). According to the calculated median risk score, 65 patients were divided into two groups (32 in the highrisk group and 33 in the low-risk group), and the clinical information is shown in Figure 3A. The PCA illustrated that patients were well divided into two clusters ( Figure 3B). The distributions of the risk score and survival time are shown in Figures 3C,D. The OS of the high-risk group was significantly worse than that of the low-risk group (p 0.0012, Figure 3E). ROC analysis of the risk model indicated that the AUC for 1, 2, and 5-year survival was 0.708, 0.815, and 0.952, respectively ( Figures 3F-H). Both of the univariate and multivariate Cox regression analyses showed that the pyroptosis-related gene signature independently predicted the prognosis of EAC patients ( Figures 3I,J).

Verification of the Gene Signature by the External Dataset
Information of 60 EAC patients from the GSE13898 dataset of GEO with available survival data was used for the validation of the gene signature. The expressions of the available differentially expressed pyroptosis-related genes are shown in Supplementary Figure S1. The patients were subdivided into the low-and high-risk groups, respectively, as described above. PCA illustrated well the separation of patients between the two groups ( Figure 4A). The distribution of the risk score and the survival time is demonstrated in Figures 4B,C. Patients in the low-risk group were with significantly higher survival rates than those in the high-risk group (p 0.003; Figure 4D). According to the ROC curve, the 1-and 2-year survival prediction models were with AUCs of 0.678 and 0.663 ( Figures 4E,F), respectively, while the 5-year survival prediction model could not be generated due to insufficient data. The risk score in our model could also serve as an Frontiers in Pharmacology | www.frontiersin.org November 2021 | Volume 12 | Article 767187 6 independent prognostic factor in the validation cohort ( Figures 4G,H).

Construction of Nomogram Based on the Gene Signature and Clinical Data
In order to more precisely predict the prognosis of EAC patients, the TNM stage was used to construct a nomogram model as shown in Figure 5A (C-index 0.764 ± 0.046). The AUCs of the nomogram for predicting 1-, 2-, and 5-year survival were 0.722, 0.884, and 1.000, respectively ( Figures 5B-D). The calibration curve indicated an ideal prediction of the nomogram ( Figure 5E). Figure 5F shows that when the nomogram-predicted probability ranged from 15% to 80%, the nomogram provided extra value relative to the treat-all-patients scheme or the treat-none scheme.

The Expressions of GSDMB and ZBP1 Within the Signature Are Downregulated by Hypermethylation
Epigenetic regulations including DNA methylation affect gene expression and modulate various cellular responses in tumorigenesis. Therefore, we further explored the mechanisms that might be involved in controlling the expressions of genes involved in the signature. We found that the RNA expressions of GSDMB (Spearman: −0.81, p 2.09e-44; Pearson: −0.80, p 5.62e-42, R 2 0.64) and ZBP1 (Spearman: −0.54, p 1.57e-15; Pearson: −0.58, p 2.59e-18, R 2 0.34) were significantly correlated with the DNA methylation status ( Figures 6A,B), whereas the association between RNA expressions of CASP1, IL1B, and PYCARD and DNA methylation was nonsignificant (Supplementary Figure  S2). Analysis by the MEXPRESS database further identified the detailed information of the methylated probes and their correlation with GSDMB ( Figure 6C) and ZBP1 ( Figure 6D) RNA expressions, suggesting that the expressions of GSDMB and ZBP1 could be regulated by epigenetic mechanisms.

Differential Expression Analysis Reveals Immune-Correlated Pathways
A total of 527 DEGs between the low-and high-risk groups were extracted according to the threshold described above. A total of 310 genes were downregulated in the high-risk group, while 217 genes were upregulated in the low-risk group. On the basis of the DEGs, GO enrichment and KEGG pathway analyses were performed. The results from GO enrichment analysis demonstrated that the DEGs were mainly associated with the regulation of cytokine production, cytokine activity, and humoral immune response pathways ( Figures 7A,B) in the TCGA cohort. KEGG pathway analysis showed that the DEGs were principally associated with the cytokine-cytokine receptor interaction and IL-17 signaling pathways ( Figures  7C,D) in the TCGA cohort. Detailed information for the deregulated pathways is shown in Supplementary Figures

S3A-D.
Similarly, the GO analysis in the GSE13898 cohort demonstrated that immune-related pathways, including neutrophil activation for immune response, neutrophilmediated immunity, and cytokine and chemokine receptor binding pathways, were deregulated between the two groups divided by the risk score (Supplementary Figures S4A, 3B). The DEGs in the GSE13898 cohort were also associated with the IL-17 signaling pathway ( Supplementary Fig. S4C, 3D).

Pyroptosis-Related Gene Signature Is Related to the Immune Status of EAC
The CMap analysis was performed to screen for small-molecular drugs that are able to revert the pyroptosis signature-related pathways, which contribute to a high-risk state. A total of 730 drugs with negative scores were identified (Supplementary Table  S3). The RAF inhibitor GDC-0879 (score: −93.76), mitogenactivated protein kinase kinase (MEK) inhibitor PD-0325901 (score: −91.95), heat shock protein (HSP) inhibitor VER-155008 (score: −89.75), mitogen-activated protein (MTOR) inhibitor torin-2 (score: −87.13), and aryl hydrocarbon receptor ligand GR-206 (score: −86.62) were the top five small-molecular drugs based on inhibition scores ( Table 1).

Pyroptosis-Related Gene Signature Is Related to the Immune Status of EAC
To explore the correlation between the selected pyroptosisrelated genes and gene-based signature with the immune microenvironment of EAC, analysis by the TIMER database for each gene was initially performed. The results indicated that ZBP1 expression was most significantly correlated with the infiltration signature of esophageal cancer, in which infiltrations of B cells (correlation coefficient 0.366, p 4.72e-07) and CD4 + T cells (correlation coefficient 0.381, p 1.41e-07) were with the most remarkable correlations ( Figure 8A, Supplementary Figure S5). In addition, somatic copy number alterations of ZBP1 were correlated with the infiltration levels of B cells, CD8 + T cells, macrophages, and dendritic cells ( Figure 8B).
The variations in the abundance of immune cell infiltration between low-and high-risk groups were further explored. The  Table  S4) and GSE13898 cohorts (Supplementary Table S5). The overview of immune cell compositions is illustrated in Figure 9A for the TCGA cohort and in Figure 9B for the GSE13898 cohort. The high-risk group in the TCGA cohort possessed significantly higher infiltration levels of M2 macrophages, activated mast cells, and eosinophils, whereas the infiltration levels of plasma cells were significantly lower ( Figures 9C,D). By contrast, the infiltration levels of memory B cells and M1 macrophages were upregulated in the high-risk group of the GSE13898 cohort, while those of naïve B cells were significantly downregulated ( Figures 9E,F).

DISCUSSION
Cell death serves as an essential barrier against the development of cancer, and pyroptosis is one of the major forms of programmed cell death (Bertheloot et al., 2021). However, the role of pyroptosis in EAC remains largely unclear. In the present study, we comprehensively evaluated the pyroptosis-related gene profiles in EAC and constructed a novel five-gene risk signature (CASP1, GSDMB, IL1B, PYCARD, and ZBP1) by LASSO Cox regression analysis. The five-gene signature showed good performance for predicting EAC prognosis in both the internal and external validation cohorts. Within the signature, GSDMB expression is distinctly correlated with the methylation status. Further enrichment analyses revealed that the DEGs between the low-and high-risk groups were correlated with immune-related pathways. The RAF inhibitor GDC-0879 and the MEK inhibitor PD-0325901 might be promising in reverting the pyroptosisrelated pathways in the high-risk EAC patients. Tumor immune microenvironment analyses indicated that high-risk patients had decreased levels of infiltrating active immune cells and higher proportions of quiescent immune-cell infiltration.
For the components within the signature, Gasdermin B (GSDMB) belongs to the GSDM family and is more broadly expressed compared to other GSDM family members (Saeki et al., 2009). The cleavage of GSDMB induced by lymphocyte-derived granzyme A triggers pyroptosis (Zhou et al., 2020). Therefore, the downregulation of GSDMB is associated with poorer prognosis of the patients. Caspase 1 encoded by CASP1 is a member of the caspase family, which is activated by inflammasomes and induces pyroptosis (Miao et al., 2011). By contrast, caspase 1 can direct T cell-independent tumor proliferation and correlates with a poorer prognosis (Zeng et al., 2018). Interleukin 1 beta (IL-1β) is a proinflammatory cytokine involved in pyroptosis. CASP-1 directly cleaves GSMD and precursor cytokines into pro-IL-1β, which initiates pyroptosis and maturation of IL-1β, respectively (Man et al., 2017). IL-1β has pro-tumorigenic effects by promoting proliferation, migration, metastasis, and angiogenesis (Gelfo et al., 2020;Rébé and Ghiringhelli, 2020). In the present study, CASP1 and IL1B upregulation is associated with a worse prognosis of EAC patients. Apoptosis-associated speck-like protein containing a CARD (ASC/PYCARD) is encoded by the PYCARD gene and contains a caspase activation and recruitment domain (CARD) for binding and facilitating the activation of caspase 1 (Bergsbaken et al., 2009). The dual role of the inflammasome adaptor PYCARD is identified in cancer cells, and therefore, PYCARD can be associated with lower cancer risks (Protti and De Monte, 2020). Z-DNA-binding protein 1 (ZBP1)-NLR Family Pyrin Domain Containing 3 (NLRP3) is critical in inducing pyroptosis by leading to cytokine maturation and GSDMD cleavage (Zheng and Kanneganti, 2020). ZBP1 expression was found to reduce tumor cell motility and chemotaxis, which decreased the potential of metastasis of tumor cells (Lapidus et al., 2007). ZBP1 stabilizes intercellular connections and focal adhesions, which suppresses breast cancer cell invasion (Gu et al., 2012). PYCARD and ZBP1 were identified as downregulated in the high-risk EAC populations. Therefore, inflammasome components might exert different effects in tumor development and progression depending on the biological context, and further investigations are needed.
Epigenetic regulation mechanisms, particularly DNA methylation, modify gene expression and regulate various cellular responses in cancer including proliferation, invasion, apoptosis, and senescence . Our study reveals that GSDMB promoter hypermethylation most notably induces decreased expression levels, indicating that methylation is essential for the regulation of pyroptosis in EAC. In recent years, epigenetic drugs are emerging, and hundreds of clinical trials are ongoing for investigating the effects of anti-DNA methylation therapies . Therefore, our results suggest that epigenetics-targeted therapy is a promising strategy for future anticancer therapeutics in part by modulating the pyroptosis-related genes. Nomograms are promising for use in clinical practice for evaluating the prognosis of EAC patients, in which the survival can be predicted using specific parameters. As indicated by the ROC curves, the nomogram demonstrates high predictive accuracy and sensitivity. Compared to the conventional TNM staging and a previously developed ferroptosis-related gene signature (AUC 0.744) in EAC (Zhu et al., 2021), the pyroptosis-related gene signature-based nomogram, which integrates gene expression profiles and clinical parameters, more effectively predicts the prognosis of EAC patients. In addition, the prognostic value of our signature is better than the DNA repair-based gene signature (AUC 0.759) in esophageal cancer . The prognostic value for 3-and 5-year survival is also higher than a recently developed signature based on nine immune-related genes for esophageal cancer (AUC 0.826) . The use of nomogram based on integrated information can facilitate the prediction of prognosis, clinical decision-making, and patient counseling (Bobdey et al., 2018).
The tumor immune microenvironment is diverge and complex, which contributes to tumorigenesis and modulates the effects of immunotherapy to a large extent. Current studies on lymphocytes in tumor immunity predominantly focus on T cells, while the protective effect of B cells has also been revealed . By contrast, mast cells have been reported to induce cancer growth (Maciel et al., 2015). Activated T cells, natural killer cells, and macrophages are potent suppressors that mediate the tumor microenvironment and exert antitumor functions (Lin et al., Frontiers in Pharmacology | www.frontiersin.org November 2021 | Volume 12 | Article 767187 2013; Nurieva et al., 2019;Cózar et al., 2021). Although some of the comparisons were not statistically different and might be contributed by the limited number of samples in both cohorts, accumulation of immune cells that promote cancer in the tumor microenvironment was generally observed in the high-risk group in both the TCGA and GEO cohorts, while the compositions of tumor-protective immune cells were reduced compared to the low-risk group. The strength of our study is that a systemic analysis was performed based on the TCGA and GEO cohorts, and the pyroptosis-related genes were assessed for the first time.
Limitations also exist in our study. Current publicly available datasets are limited in both number and size, and therefore, validation of our prediction model in large-scale EAC cohorts could be performed in future studies. In addition, based on the information of our study, further in vitro and in vivo studies could be conducted to evaluate the function and mechanisms of pyroptotic regulation in EAC. Despite the limitations, our study has provided a comprehensive overview of pyroptosisrelated gene profiles in EAC.
In summary, we identified differentially expressed pyroptosisrelated genes and developed a novel five-gene pyroptosis signature that significantly correlates with the survival of EAC patients. The pyroptosis-based signature is an independent prognostic factor and performs better than the TNM stage, which is promising for clinical application. Moreover, GSDMB expression is notably correlated with methylation status, and the signature is related to antitumor immunity in the tumor microenvironment. Modulating pyroptosis, epigenetic mechanisms, and immune microenvironment by drug discovery might be promising for improving the prognosis of patients. Further studies exploring the regulating patterns are warranted.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
RZ and SH acquired the data, performed the analysis, and wrote the manuscript. XQ was involved in the analysis, data interpretation, and revision of the manuscript. ZZ and HW participated in data analysis. HC, WS, and LJ were involved in study design, supervision, and acquiring funding.