Single-Cell RNA Sequencing Reveals the Immunological Profiles of Renal Allograft Rejection in Mice

Allograft rejection is a common immunological feature in renal transplantation and is associated with reduced graft survival. A mouse renal allograft rejection model was induced and single-cell RNA sequencing (scRNA-seq) data of CD45+ leukocytes in kidney allografts on days 7 (D7) and 15 (D15) after operation was analyzed to reveal a full immunological profiling. We identified 20 immune cell types among 10,921 leukocytes. Macrophages and CD8+ T cells constituted the main populations on both timepoints. In the process from acute rejection (AR) towards chronic rejection (CR), the proportion of proliferating and naïve CD8+ T cells dropped significantly. Both B cells and neutrophils decreased by about 3 folds. On the contrary, the proportion of macrophages and dendritic cells (DCs) increased significantly, especially by about a 4.5-fold increase in Ly6cloMrc1+ macrophages and 2.6 folds increase in Ly6cloEar2+ macrophages. Moreover, myeloid cells harbored the richest ligand and receptor (LR) pairs with other cells, particularly for chemokine ligands such as Cxcl9, Cxcl10, Cxcl16 and Yars. However, macrophages with weak response to interferon gamma (IFNg) contributed to rejection chronicization. To conclude, reduction in CD8 T cells, B cells and neutrophils while increasing in Ly6cloMrc1+ macrophages and Ly6cloEar2+ macrophages, may contribute significantly to the progress from AR towards CR.


INTRODUCTION
Rejection of the transplanted kidney in humans is still a major cause of morbidity and mortality. Despite current immunosuppressive therapy, AR develops in about 10-12% of transplant patients (1). Besides, the long-term survival rate is not satisfactory, as nearly half of the transplanted kidney function gradually loses within 10 years (2). scRNA-seq has enabled the profiling of specific cell populations at the single-cell level, offering an unprecedented opportunity to define cell types and states comprehensively with molecular precision. However, due to unpredictability of tissue availability, limited size of kidney biopsy samples, challenges are raised in wide application of scRNA-seq to dissect immunological network and dynamic changes involved in human kidney allograft rejection. The mouse genome has been well characterized and the proteincoding regions of mouse and human genomes are 85% identical (3). Thus, mouse models of diseases are available for scRNA-seq to reveal the characteristics of immunity (4)(5)(6). Of them, the fully mismatched mouse model of renal transplantation is a better model that closely relates to that of human both technically and pathologically (7). Studying the evolution of the pathology of rejection in mouse kidney allografts, interstitial infiltrate, major histocompatibility complex (MHC) induction, and venulitis peaked at day 7-10, then stabilized or regressed by day 21 (8). By scRNA-seq, immune cells isolated from allografts undergoing AR on D7 to CR on D15 following kidney transplantation should gain new insights into the immunological profiling, dynamic changes from AR to CR during renal transplant rejection.

Mouse Kidney Transplantation
Ectopic renal transplantation was performed by transplanting kidneys from male BALB/c mice (H-2d; 8 weeks old) into male C57Bl6/J recipients (H-2b; 8 weeks old) as previously described (9,10). Mice were euthanized on D7 after the operation. Animal studies were performed according to the Department of Health (Hong Kong) guidelines in the Care and Use of Animals and the experimental protocol was approved by the Animal Experimentation Ethics Committee at the Chinese University of Hong Kong.

Histology
Kidney tissues were fixed in formalin and embedded in paraffin, sectioned into 3 µm slices for hemotoxylin & eosin (HE), periodic acid-Schiff (PAS) and Masson staining.

Serum Enzyme-Linked ImmunoSorbent Assay (ELISA)
Blood was collected in serum separator tube by cardiac puncture from WT C57Bl6/J mice, kidney allograft recipients sacrificed on D7 and D15 after kidney transplantation. Samples were allowed to clot for 30 min at room temperature before centrifugation for 15 min at 1,000g, 4°C. Remove serum, aliquot, and store samples at −80°C. Serum IgG and IgM ELISA (Cat# ab157719, ab215085; Abcam, Quantikine, USA) was performed as per the manufacturers' recommendations.

Fluorescence-Activated Cell Sorting (FACS)
After sacrifice, mice were flashed with 0.9% saline through heart to remove the peripheral blood remained in the kidney allografts. Kidney graft tissues were then collected and digested using Blenzyme 4 (Roche Inc., Indianapolis, IN). Kidney suspensions were pooled and stained with Live/dead Aqua Fluorescent Reactive Dye (Life Technologies) and PE conjugated antimouse CD45 antibody (eBioscience). DAPI -CD45 + cells were FACS-sorted on (BD Influx ™ cell sorter) with a 100 µm nozzle and collected ( Figure S1A).

A 10× Sample Processing, cDNA Library Preparation and Single Cell Mapping
Cells were assessed for viability and counted using Countess II FL Automated Cell Counter, loaded into individual wells of 10× Chromium Single Cell chip. Single cells were then encapsulated into Gel in-beads emulsions by 10× Chromium Controller. Chromium ™ Single Cell 5' Reagent Kits (v1 chemistry) were used to perform reverse transcription on Gel in-beads emulsions, followed by cDNA clean up and amplification. The doublestranded cDNA was then gone through enzymatic fragmentation, adapter ligation, index PCR and SPRI select size selection per manufacturer's protocol. Library size and concentrations were determined by Qubit, quantitative PCR and Bioanalyzer assays. scRNA-seq data were mapped to the mouse (mm10) reference genome and quantified using CellRanger 2.2.0 supplied by 10×.

D15 CD45 + Leukocytes Extraction and Data Sets Integration
D15 data sets, GSM4761000 and GSM4761003, contained in dataset GSE157292 were downloaded from the Gene Expression Omnibus (GEO). Cells from both D7 and D15 datasets were included with the exclusion criteria (1): cells with greater than 6% mitochondrial genes (2), cells expressing less than 200 genes or greater than 5,000 genes, and (3) cells with less than 1,000 unique molecular identifier (UMI) counts or greater than 40,000 UMI counts. All these cells were removed to filter out low-quality cells and doublets. Next, Uniform Manifold Approximation and Projection (UMAP) were used for dimension reduction and clustering by Seurat 3.0 version on dataset GSE157292. Cell clusters expressing gene Ptprc (cluster 0: 7, 9, 11, 12, 15:19) were then extracted as D15 dataset (Supplementary Figure 2). Normalization, feature selection, scaling and principal component analysis were then performed separately for D7 and D15 datasets before integration. Again, UMAP were used for dimension reduction and clustering.

Single Cell Marker Identification and Differential Expression
Marker genes were determined for each cluster by Wilcoxon Rank Sum test within the FindAllMarkers function using genes expressed in a minimum of 10% of cells and average natural log-transformed fold change threshold of 0.25. We defined a gene to be a cluster-specific marker if it was differentially upregulated in a cluster as compared to the remaining clusters. To identify differential expression genes (DEGs), we used Seurat FindMarkers function using genes expressed in a minimum of 10% of cells and with an adjusted p value (adj_p value) threshold of 0.05.
Gene Ontology (GO) and Gene Set Variation Analysis (GSVA) GO enrichment analysis was performed on Metascape (metascape.org), gene-set enrichment score was calculated with GSVA software from Bioconductor (version 1.28). Marker genes or DEGs with average log 2 fold change (avg_log 2 FC) >0.59, adj_p value <0.05 were used for GO. Gene-sets were downloaded from the Mouse Genome Informatics (http://www.informatics.jax. org/). Raw counts file was used as input data. R/Bioconductor Limma package was used to test whether the gene sets were either upregulated on D7 or D15.

Single-Cell Regulatory Network Inference and Clustering (SCENIC)
SCENIC was performed following pySCENIC (11). Activity score of each regulon was calculated in each single cell through the area under the recovery curve to evaluate whether a regulon was enriched at the top ranking of cells. Transcription factors (TFs) motif enrichment analysis was then performed to identify the direct targets (regulons) and score the activity of the regulons (AUCell score). AUCell score was then projected onto the UMAP or scaled to plot heatmaps.

Pseudotime Analysis
The original gene expression matrix of mononuclear/ macrophage cell populations (C1, C2, C6, C7, C12 and C16) were used as the input files for creating the Monocle 2.0 CellDataSet. The DDRTr ee m et ho d w a s u s e d fo r dimensionality reduction. Batch effect was removed by the residualModelFormulaStr function.

Ligand-Receptor (LR) Interaction
We conducted cell-cell interaction analysis utilizing cellphonedb function curated by the CellPhoneDB database (12). The significant cell-cell interactions were selected with p value <0.01. Heatmaps of significant interaction pairs were plotted using mean values of interacting partners. Mean value refers to the total mean of the individual partner average expression values in the corresponding interacting pairs of cell types.

Statistics
Data were expressed as the mean ± SD and analyzed using student t test and one-way analysis of variance, followed by Newman-Keuls multiple comparison test from GraphPad Prism 6.07 (GraphPad Software Inc, San Diego, CA, USA). p < 0.05 was considered statistically significant.

AR Progressed Towards CR From D7 to D15 During Mouse Renal Allograft Rejection
After tissue harvesting, histological features of kidney allografts from WT and kidney allograft recipient mice were examined by HE, PAS and Masson staining ( Figure 1A). Native kidneys do not exhibit tubular injury. Diffuse mononuclear cell infiltrates throughout the renal parenchyma, severe arteritis, necrotic tubules and tubulitis were observed in allografted kidneys on both D7 and D15 following transplantation. However, Masson staining also detected that moderate fibrosis was developed on D15 after the operation, whereas almost no collagen fibers were observed on D7 ( Figure 1A). Thus, histologic evidence indicated that it was AR on D7 and a mix rejection pattern on D15. Interestingly, although serum levels of IgM were markedly elevated on both 7 and 15 days after kidney transplantation, an increase in serum levels of IgG was only observed on D15 ( Figures 1B, C), suggesting a switch from the early acute T cell-mediated rejection to the late mixed T cells and antibodymediated rejection.
To observe polarization of macrophage subtypes, expression of M1 and M2 markers were plotted in Supplementary Figure 4. Cluster 1 and 16 highly expressed M2 markers Mrc1, Clec10a and Selenop. Besides Mrc1, they also expressed another scavenger receptor Stab1 ( Figure 3D), suggesting a role in scavenging debris/excess exracellular matrix. Ly6C hi macrophages (cluster 6) significantly expressed M2 marker Chil3; Ly6C inter macrophages (cluster 2) significantly expressed M1 markers Il1b and Nos2; Ly6C lo macrophages (Cluster 7) also highly expressed Il1b; Tnf was widely expressed by all macrophage clusters, except cluster 12, a group of cells that do not express any of the above-mentioned polarization markers. Then, marker genes in each macrophage cluster were used for GO enrichment analysis (Supplementary were selected to compare the enrichment scores among macrophage subtypes by gene set variation analysis (GSVA, Figure 3E). In general, clusters 1, 2, 6, 7 contributed to inflammatory response during rejection, whereas cluster 12 seemed to be quite quiescent with no obvious functions. Cluster 16 is a group of proliferatin macrophages. Functionally, Ly6C lo Mrc1 + macrophages (cluster 1) may contribute to chemokine signaling pathway, receptor-mediated endocytosis, activation of ERK1 and ERK2 cascade. Ly6c inter Nos2 + macrophages (cluster 2) significantly correlated with response to IFNg, antigen processing and presentation, reactive oxygen species metabolic pathway, thus were annotated as IFNg induced cells (IFNIC). Ly6c hi Chil3 + macrophages (cluster 6) were enriched in pathways response to interferon beta and il1b, chemokine/cytokine signaling pathway, phagocytosis, and wound healing. Ly6C lo Ear2 + macrophages (Cluster 7) showed activation of I-kappaB kinase/NF-kappaB signaling pathway. We next applied SCENIC analysis to assess which transcription factors (TFs) were responsible for the different gene profiles and functions among macrophage clusters (Supplementary Figures 6A-H) (11). SCENIC identified Maf, Tcfl2, Mitf and Nr1h3 as main candidate TFs underlying the specific gene expression in cluster 1 Ly6C lo Mrc1 + macrophages. Usf1, Spi1 and M1 polarization associated TFs Hif1a, stat1 were activated in IFNIC (cluster 2), whereas Cebpa, Xbp, Irf5 and Jund were activated in Ly6C hi Chil3 + macrophages (cluster 6). Cebpb, Etv6, Rel, Elf4, Elf2, Creb1, Elk3 and Smarcc2 were underlying the specific gene expression in Ly6C lo Ear2 + macrophages (cluster 7). TFs related to proliferation, including E2f2, E2f7, E2f8 and Mybl2, were specifically activated in cluster 16.

Blunted Response to Interferon Gamma Contributed to AR Chronicization
To observe the functional changes of macrophages from AR towards CR, differential expression genes (DEGs) comparing total macrophages (cluster 1, 2, 6, 7, 12, 16) on D7 versus D15 were calculated. 292 DEGs with avg_ log 2 FC >0.59, adj_p value <0.01 were upregulated on D7 whereas 197 were upregulated on D15 (Supplementary Figure 7A). GO enrichment analysis results were shown in Figure 5A. In the meanwhile, GSVA analysis was conducted (Supplementary Figure 7B). Both indicated that macrophages infiltrated in kidney allografts on D7 following transplantation were more activated in the pathway response to IFNg than those on D15 ( Figures 5B, D; Supplementary Figure 7B). Representative DEGs enriched in response to IFNg (Nos2, Socs1, Stat1, Cxcl10, Mif, etcetera) were marked in Figure 5D. Interestingly, we found a significant number of ribosomal protein genes that  were highly and consistently expressed on D15 than D7 ( Figures 5C, E), resembling a tolerized gene signature reported by others (13).

NK Cells
A small population of renal allograft NK cells were identified by NK-specific markers ( Figures 2B, 8A), expressing Eomes, a T-box transcription factor expressed by activated CD8 T cells and NK cells ( Figure 8B) (19). They also significantly expressed mRNA for proinflammatory cytokine Ifng and chemokine Ccl3 (MIP-1a), Ccl4 (MIP-1b), Ccl5, Xcl1, which attract effector lymphocytes and myeloid cells toward inflamed tissues ( Figure 8B) (20). Besides, NK cells expressed cytotoxic marker genes Prf1, Gzma, Fasl, playing roles in the regulation of natural killer cell mediated cytotoxicity and positive regulation of cytokine production ( Figure 8C). SCENIC identified Irf8 as main candidate TF underlying the specific gene expression in NK cells ( Figure 6D, cluster 13), which regulated natural killer cell development (21).

B Cells
We detected 469 B cells that formed two clusters ( Figure 9A), with the proportion dropping from 7.85% on D7 to 1.38% on D15 ( Figure 2C, Supplementary Table 1). Cluster 10 was the major B cell population at both peak and regression phases, characterized by Cd79a, Cd79b and Siglecg ( Figure 2B). Cluster 19 (plasma B cell) was marked by Jchain and Eaf2. Besides, cluster 10 expressed MHCII molecules and Ighd at higher levels, whereas cluster 19 (plasma B cells) highly expressed genes encoding immunoglobulins, including Ighg2b, Igha, Iglc1, Iglc2, Iglv2, and Ighj3, etcetera ( Figure 9B). The functions of B cell cluster 10 included translation, antigen processing and presentation of peptide antigen via MHC class II, response to interferon-gamma. The functions of plasma cells (cluster 19) included response to endoplasmic reticulum stress and immunoglobulin production ( Figure 9C). SCENIC identified Foxo1, Runx1, Jun and Tgif1 as main candidate TFs underlying the specific gene expression in B cell cluster 10, Pbx1 and Xbp1 as main candidate TFs underlying the specific gene expression in plasma cell cluster 19 ( Figure 9D).

Granulocytes
Neutrophils are usually the first leukocytes to infiltrate into transplanted organs and are a well-established marker of transplant injury (22). Only one neutrophil cluster (cluster 15) with 228 cells was identified in our scRNA-seq data, constituting <5% of the immune spectrum ( Figure 10A, Supplementary Table 1). The proportion of neutrophils dropped about 3.5 folds from D7 to D15 (Supplementary Table 1). Featured higher expression of mRNA for pro-repair genes Mmp8, Mmp9 and Arg2 were detected in neutrophil cluster, resembling a cluster of pro-repair subset reported by others ( Figure 10B) (6). Besides, Il1b was expressed at particular high levels in neutrophils, contributing to the pathogenesis of immune response during inflammation. Basophils were also present in the majority of renal biopsies showing T cellmediated rejection and were absent in biopsies without rejection (23). A small population of basophils (cluster 18, <1%) was identified by expressing basophil marker Mcpt8, Gata2 ( Figure 2B) and activation markers, including type 2associated cytokines Il4, Il13, and proinflammatory marker Il6 ( Figure 10B, Supplementary Table 1). CSF1, a critical growth factor for macrophage development, was predominantly expressed by granulocytes, including neutrophils and basophils, among all the immune cells ( Figure 10C). SCENIC identified Stat1, Ets2 and Cebpb as main candidate TFs underlying the specific gene expression in neutrophils, Gata1, Gata2 and Fosl1 as main candidate TFs underlying the specific gene expression in basophils ( Figure 10D).

Complex Intercellular Communication Networks Among Immune Cells
Next, we used CellPhoneDB to identify LR pairs and molecular interactions among the major cell types (12,24). (Supplementary Figures 10A, B), followed by macrophages ( Figures 12A-F). Notably, both macrophages and cDCs commonly expressed relatively high levels of chemokines (e.g., Cxcl9, Cxcl16, Yars and Cxcl10), whereas, the corresponding receptors were widely expressed in other immune cells, suggesting that these chemokines play significant roles in interactions within myeloid cells or between myeloid cells and other cells ( Figures 12A-F, Supplementary Figures 11A, B). Compared to other macrophage clusters, Ly6c lo Ear2 + macrophages exhibited significantly less and low means of LR pairs, lacking expression of receptors including Ccr1, Ccr2, Ccr5 and Ccr8 ( Figure 12D). Ly6c lo Mrc1 + macrophages (cluster 1, 16) and IFNIC (cluster 2) specifically expressed ligand Axl and interacted with other cells through Axl-Il15ra, Axl-Pros1 pairs ( Figures 12A, C, F, Supplementary Figure 11C). It was reported that Axl augmented intragraft differentiation of certain types of proinflammatory macrophages and enhanced early allograft inflammation (13).

DISCUSSION
Here, we presented, for the first time to our knowledge, a singlecell resolution catalog of immune cells profiling from AR towards CR during mouse renal allograft rejection. It is reasonable to integrate our single-cell data of D7 with D15 data extracted from GSE157292, as kidney transplant models were both built by transplanting kidneys from BALB/c into B6 mice. Though Dangi et al. bilaterally nephrectomized the native kidneys while we left them in situ, the histological outcomes on D7 are comparable between two mouse models (8). In addition, two scRNA-seq studies were both performed on 10× genomics platform. There were several novel findings from this study. First, we found a significant change in the proportion of immune cells in the grafted kidneys between AR and mix rejection. Generally, the proportion of macrophages increased from 21.11 to about 50%, including Ly6C lo Mrc1 + macrophages (clusters 1, 16) and Ly6C lo Ear2 + macrophages (cluster 7). Pseudotime analysis indicated that both Ly6C lo Mrc1 + and Ly6C lo Ear2 + macrophages were at the terminal differentiation stage. To the opposite, CD8 T cells dropped from 51.44 to about 28%, mainly due to the decrease in proliferating CD8 T cells (cluster 8) and naïve CD8 T cells (cluster 4). In addition, B cells and neutrophils also showed significant decrease. Second, the macrophage polarization with distinctive functions may contribute to renal allograft progression or regression. Ly6C hi monocytes recruitment could be induced by type I IFN signaling (25,26). During acute renal ischemic-reperfusion injury, it was found that the CD11b + Ly6C hi macrophage population was associated with the onset of renal injury and increased expression of proinflammatory cytokines, whereas the Ly6C inter population had a distinct phenotype of wound healing and the Ly6C lo population exhibited a profibrotic phenotype as previously reported (27). These findings indicate that the macrophage phenotypes may contribute significantly to the progression or regression of renal allograft outcomes. Interestingly, by GO enrichment and GSVA analysis, we also showed that Ly6C hi Chil3 + (coding Ym1) macrophages (cluster 6) responding to type I IFN and Il1b were enriched in wound healing and phagocytosis, contributing to the resolution of inflammation and tissue repair. Ly6C lo macrophages could be further clustered into Ly6C lo Ear2 + macrophages (Cluster 7) with weak chemoattraction and Ly6C lo Mrc1 + macrophages (Clusters 1 and 16) acquiring a resident macrophage-like gene signature, both increased in the process towards CR. IFNIC (cluster 2), induced by IFNg and ROS, also gained high GSVA score in antigen processing and presentation. Pseudotime trajectories indicated that IFNIC may be a heterogeneous subtype. These findings suggested that although Ly6C is used to classify macrophages, functions of macrophage subtypes could be altered during the renal allograft progression or regression. This can also be observed in other disease models. For example, in unilateral ureteral obstruction model, kidney infiltrating immature pro-inflammatory Ly6C hi macrophages can progressively differentiate into mature profibrotic or M2 Ly6C lo macrophages (28,29). In mouse cardiac transplantation, recipient derived suppressive Ly6C lo macrophages show inhibitory effect on CD8 + T cell immunity while promoting Treg cell expansion. Besides, immunogenic Ly6C hi macrophages could develop into tolerogenic Ly6C lo macrophages under the action of neutrophil derived CSF1 (30,31). Despite the generally pro-inflammatory effects of Ly6C hi monocytes/macrophages, it has also been reported that a cluster of bone marrow originated Ym1 + Ly6C hi monocytes exhibited immunoregulatory and tissuereparative phenotypes during the recovery phase of colitis (32). In addition, Ly6C lo macrophages may also be pro-inflammatory and induce kidney injury in acute kidney injury mice (33). Thus, macrophages are highly heterogeneous and participate in the disease inflammation or resolution. Third, the present study proved IFNg as an important inducer of macrophage and cDC activation during mouse renal allograft rejection. T lymphocytes and NK cells significantly expressed IFNg ( Figure 2B). Macrophage subtypes and cDC clusters were also enriched in pathway response to IFNg (Supplementary Figures 6A-H). Importantly, both DEGs enrichment and GSVA analysis also showed that IFNg response decreased from the inflammation peak on D7 to D15. In addition, we found a wide range of LR pairs among immune cells. Granulocytes and macrophages interacted most extensively with other cell types. Ligation by IFNg induced chemokines like Cxcl9, Cxcl16 and Cxcl10 in macrophages and cDC mediated the activation of other cells. Indeed, IFNg is needed for transplant rejection by inducing MHC products and pro-inflammatory molecules in both human and mouse kidneys (34)(35)(36)(37)(38)(39). In vitro study also proves that IFNg can stimulate the maturation of DC by upregulating MHCII molecules and costimulatory molecules CD40 and CD80 (40). Thus, decrease in IFNg response by renal macrophages and cDCs may contribute to inflammation regression during renal allograft regrection.
So far, only limited information on the use of scRNA-seq in murine kidney transplantation has been published. Anil Dangi et al. reported whole transcriptomics of mouse kidney transplants on day 15 in a single cell level and identified 13 immune cell clusters (13). In the present study, we developed further in-depth analysis of these immune cell types with comprehensive classification and description of cell subtypes including their specific marker genes, functional characteristics, transcriptional regulation, and the proportional dynamics from the renal allograft rejection peak on D7 to the regression phase on D15. In naïve mouse kidneys, monocytes/macrophages, neutrophils, DCs, B and T lymphocytes, NK cells could all be observed (16,41,42). Among macrophage subtypes, infiltrating Ly6C hi Chil3 + macrophage and Ly6C lo patrolling macrophages are commonly found in normal kidneys (16,42), whereas IFNIC and Mrc1 + macrophages are more common to disease models like unilateral ureteric obstruction model (42). Study of kidney allograft biopsies from human recipients undergoing AR and CR also identified immune populations like those in mice, including monocytes, B cells, T cells, plasma cells (43)(44)(45), therefore independently supporting a role of these cells in kidney rejection in humans. However, due to the limited size of human biopsy samples and the large percentage of parenchymal cells, further analysis of immune cell subtypes is lacking. Thus, in-depth analysis of immune cells in human kidney biopsies is needed. In summary, results from this study may provide more comprehensive insight for our better understanding of the immune cell profiling from AR towards CR and may also provide valuable information for the future studies in immunological targeting and tolerance.

DATA AVAILABILITY STATEMENT
Raw scRNA-seq data (fastq files) and processed unique molecular identifiers count matrices used in the current study are publicly available on NCBI GEO (GSE166775).

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Experimentation Ethics Committee at the Chinese University of Hong Kong.

AUTHOR CONTRIBUTIONS
QS, JHC, HJ, and HL conceptualized the study. QS, YW, JYC, and XH performed experiments. QS and LM analyzed bioinformatics. JYC, HJ, HL, and ST provided reagents and animals. QS drafted the manuscript. JYC, HJ, and HL edited and finalized the paper. All authors contributed to the article and approved the submitted version.