Single-cell RNA-seq unveils critical regulators of human FOXP3 + regulatory T cell stability

The heterogeneity and plasticity of T lymphocytes is critical for determining immune response outcomes. Functional regulatory T (Treg) cells are commonly characterized by stable FOXP3 expression and have reported to exhibit heterogeneous phenotypes under inﬂammatory conditions. However, the interplay between inﬂammation and Treg cell suppressive activity still remains elusive. Here, we utilized single-cell RNA sequencing to investigate how human Treg cells respond to the pro-inﬂammatory cytokine interleukin-6 (IL-6). We observed that Treg cells divided into two subpopulations after IL-6 stimulation. TIGIT (cid:1) unstable Treg cells lost FOXP3 expression and gained an effector-like T cell phenotype, whereas TIGIT + Treg cells retained robust suppressive function. Single cell transcriptome analysis revealed a spec-trum of cellular states of IL-6-stimulated Treg cells and how cytochrome P450 family 1 subfamily A mem- ber 1 (CYP1A1) is a crucial regulator of Treg cell suppressive capability and stability. CYP1A1-deﬁcient human Treg cells developed a Th17-like phenotype after IL-6 stimulation. Our ﬁndings implicate CYP1A1 as a previously unidentiﬁed regulator of Treg cells that may have target potential for clinical application for biotherapies. (cid:1) 2020 Science China Press. Published by Elsevier B.V. and Science China Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


a b s t r a c t
The heterogeneity and plasticity of T lymphocytes is critical for determining immune response outcomes. Functional regulatory T (Treg) cells are commonly characterized by stable FOXP3 expression and have reported to exhibit heterogeneous phenotypes under inflammatory conditions. However, the interplay between inflammation and Treg cell suppressive activity still remains elusive. Here, we utilized singlecell RNA sequencing to investigate how human Treg cells respond to the pro-inflammatory cytokine interleukin-6 (IL-6). We observed that Treg cells divided into two subpopulations after IL-6 stimulation. TIGIT À unstable Treg cells lost FOXP3 expression and gained an effector-like T cell phenotype, whereas TIGIT + Treg cells retained robust suppressive function. Single cell transcriptome analysis revealed a spectrum of cellular states of IL-6-stimulated Treg cells and how cytochrome P450 family 1 subfamily A member 1 (CYP1A1) is a crucial regulator of Treg cell suppressive capability and stability. CYP1A1-deficient human Treg cells developed a Th17-like phenotype after IL-6 stimulation. Our findings implicate CYP1A1 as a previously unidentified regulator of Treg cells that may have target potential for clinical application for biotherapies.
Ó 2020 Science China Press. Published by Elsevier B.V. and Science China Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

Introduction
CD4 + CD25 + regulatory T (Treg) cells are essential for immune homeostasis. Transcription factor forkhead-box-P3 (FOXP3) is specific expressed in Treg cells and play a crucial role in maintaining the expression of Treg immunosuppressive signature genes [1]. FOXP3 controls the expression of a range of immune modulatory genes, which distinguishes Treg cells from CD4 + FOXP3 À conventional T cells [2,3]. In addition, Treg cells exhibit both instability and plasticity upon cytokine and antigen challenge [4][5][6]. It has been depicted by recent studies that immunosuppressive capability of Treg cells can be depressed under specific environmental signals while inflammatory effector gene expression is induced [7][8][9][10][11].
It has become clear that extensive heterogeneity exists in Treg cells with specific functions [12,13]. TIGIT as a co-inhibitory molecule, is highly expressed in Treg cells and directly inhibits T cell activation and proliferation. TIGIT high expressed Treg cells may define an activated functionally Treg cell subset [14]. As a pleiotropic cytokine, IL-6 has various biological functions, including hematopoiesis, inflammation and oncogenesis [15]. Elevated IL-6 production in pathogenic and inflammatory tissue has been suggested to reduce the suppressive capacity of Treg cells and promote the differentiation of IL-17-producing Th17 cells [16][17][18]. The heterogeneity and plasticity of Treg cells become more and more clear [19][20][21]. But the relationship between inflammation and Treg heterogeneity are still poorly understood. Therefore, it is our interest to unveil the heterogeneity within Treg populations and to understand how Treg subsets respond to inflammatory cytokine stimuli at the single cell level [13,22].
Single-cell analysis of lymphocytes opens the way for deep understanding of these cells in different tissues. Through single cell sequencing, it has been recently reported that TCR signaling can shape Treg differentiation and effector function [23]. Single-cell transcriptome analysis has been used to several research areas including tumor immunity, immune cell development, dendritic cell bimodality and hematopoiesis [24][25][26][27][28]. Single-cell transcriptome analysis can also allow for further insight into the potential regulators of immune cell function. For instance, single-cell RNA sequencing experiments have shown that Th2 cells can inhibit effector T cell function and proliferation by Cyp11a1 mediated steroidogenesis [29]. Single-cell RNA-seq in Th17 cells unveiled four novel functional genes including Gpr65, Plzp, Toso, and Cd5l, and verified their substantial effects on autoimmune diseases [30]. In tumor-infiltrating T cells, LAYN was identified as a negative regulator of antitumor immunity, where elevated expression was found associated with poor prognosis in hepatocellular carcinoma patients [26].
Here, we defined the phenotype and functional heterogeneity of Treg populations after treatment with the proinflammatory cytokine IL-6 [15,31]. We identified two subpopulations of Treg cells after IL-6 stimulation. An unstable Treg subpopulation had reduced expression of FOXP3 and Treg associated markers, elevated TNF-a and IL17A expression, and increased glycolysis and proliferation. A rather stable Treg subpopulation showed higher Treg-associated marker expression. We propose that, the immunosuppressive function of partially stable Treg cells can be disarmed allowing for adoption of effector functions; stable Treg cells are sensitive to inflammatory signals and infiltrate into inflamed tissue to prevent excessive inflammation in the presence of IL-6. Moreover, we found and validated CYP1A1 as a new regulator of Treg cell function and stability. Knockdown of CYP1A1 significantly impaired Treg cells immunosuppressive capability of Treg cells and increased the expression of inflammatory cytokines. Therefore, CYP1A1 could be a novel therapeutic target to regulate Treg cells function under inflammation.

Ethical approval
This study was approved by Institutional Review Board on Bioethics and Biosafety of BGI (Research Ethics Committee reference number FT 15237-2). All participants provided written informed consent according to regulation of the Declaration of Helsinki.

Single-cell separation, library construction and RNA-seq
We collected in vitro cultured human Tregs with suspension buffer using FACS by cmFDA at the concentration of 10 cells per lL, and then prepared mRNA SMART-seq libraries using the microwell full-length mRNA amplification and library construction system (MIRALCS). We selected 96 amplified cDNAs according to the cutoff cycle threshold (Ct) < 15, melting temperature (Tm) > 88, and transferred the target cDNAs into 96-well plates. Tn5 transposase-based library construction was processed manually in a 96-well plate. The amplified cDNAs and libraries were quantified. Finally, samples were sequenced using the BGISEQ-500 system (Complete Genomics TM sequencing technologies).

Sample filtering
Sample and gene filter pipeline were performed on all single cells using in-house scripts. To ensure whether the sequencing depths were sufficient, we carried out saturation analyses on 10 random selected samples and found that all samples reached saturation in gene detection with about 3 million raw reads used when we chose an FPKM cutoff at 1. Single cells with raw reads < 3 million were removed. After filtration the mean clean reads of all effective cells was up to 20.6 million. Moreover, quality control was included: (1) aligned reads rate > 70%, (2) percentage of ERCC spike-in mRNA reads < 10%, (3) Pearson correlation coefficient between ERCC spike-in mRNAs counts and ERCC spike-in mRNA expression > 0.6, (4) no 5' to 3' Bias, (5) no GC content Bias, (6) detected gene number in individual cells > 5000 (FPKM > 1). Following these criteria, 82 effective cells were kept to carry out batch correction and further gene expression analysis.

Gene filtering and batch correction
We generated two gene expression matrices of protein coding genes based on different gene filtering criterion. The strict filtering gene expression matrix required gene expression FPKM beyond 8 in at least 20% of the cells, and 3885 genes were retained. Likewise, the loosely filtering gene expression matrix required gene expression FPKM beyond 1 in at least 5% of the cells, and after filtration 9683 genes were retained. As all 82 cells were processed in the same chip, and the library was constructed in the same 96-well plate, batch correction was performed on these two filtered gene matrices using the COMBAT package.

Weighted gene co-expression network analysis (WGCNA)
According to genes expressed at more than 20% and FPKM more than 5, in total 4795 genes were selected for WGCNA. The gene coexpression networks were constructed by WGCNA R package with b = 8 and a minimum module size of 25 genes.

Pseudotime analysis
Pseudotime analysis can order different generations of Treg cells through biological processes based on transcriptional similarities. Embeddr R package was applied to pseudotime analysis with default settings [36].

Gene set enrichment analysis (GSEA)
Unbiased analysis of gene signature enrichment analysis was made by javaGSEA application with default settings [37].

Flow cytometry
All cells were washed in PBS with 2% FBS, then stained with antibodies and fixable viability dye eFluor 780 Intracellular staining of FOXP3 was performed according to the manufacturer's instructions (eBioscience). For cytokine detection, Treg cells were stimulated with PMA (50 ng/mL), ionomycin (1 mmol/L), Golgi Stop and Golgi Plug for 4 h before cytokine detection. All samples were processed on a Fortessa flow cytometer (BD Biosciences) and data were analyzed by FlowJo software (TreeStar).

RNA preparation and immunoblotting
Total RNA extraction and qRT-PCR analysis were carry out according to the manufacturer's specification (SYBR Green; TaKaRa). For immunoblotting, Treg cells were lysed on ice for 30 min in RIPA buffer containing protease inhibitor (1:100, P8340; Sigma-Aldrich), 1 mmol/L NaF, and 1 mmol/L PMSF. Loading buffer was added to supernatants and boiled at 100°C for 10 min.

In vitro suppression assay and Treg cell proliferation assay
Conventional T cells (CD 4+ CD 25-CD 127+ ) were sorted and stained with CellTrace Violet (Invitrogen). Labeled conventional T cells were co-cultured with Treg cell for 3 days. The proliferation of CYP1A1-depleted Treg cells was also analyzed using the same protocol. Flow cytometry was performed to detect conventional T cells proliferation. Treg cells isolated from cord blood were also activated with anti-CD3/CD28 antibodies and treated with CYP1A1 inhibitor Alizarin (6.2 lmol/L, Selleck) for 48 h.

Statistical analysis
Statistical significance was determined using GraphPad Prism with Student's t test as specified in figure legends. Results are expressed as means ± s.d. *P < 0.05, **P < 0.01, ***P < 0.001, NS, not significant, as determined by two-tailed, unpaired Student's t-test.

FOXP3 expression is heterogeneous after IL-6 treatment
IL-6 is a known proinflammatory cytokine that can inhibit the immunosuppressive capability of Treg cells. To examine different responses of Treg cells to IL-6 exposure, we sorted naïve Treg cells (CD4 + CD25 + CD127 dim CD45RA + CCR7 + ) from human cord blood (Fig. 1a) followed by stimulation for 3 days with anti-CD3/CD28 antibody, IL-2 plus IL-6. According to FOXP3 expression, Treg cells were characterized into unstable and stable Treg cells (Fig. 1b). Unstable Treg cells showed significant decrease FOXP3 expression compared to their stable counterparts. We also sorted Treg cells (CD4 + CD25 + CD127 dim ) from adult peripheral blood followed by stimulation for 3 days with anti-CD3/CD28 antibody, IL-2 plus IL-6. We found Treg cells heterogeneity with similarity from Treg derived from cord blood (Fig. S1 online). This observation indicated gene expression difference between stable and unstable Treg cells.

Single-cell RNA sequencing reveals Treg cell heterogeneity
Since there is only a limited number of a surface marker to distinguish between unstable and stable Treg cells, we performed single-cell RNA-seq to dissect Treg cell heterogeneity. We performed single-cell full-length transcript amplification for sequencing by MIRALCS (Fig. S2a online) [38]. Then we performed rigorous quality control to remove cells with poor gene coverage and expression (Fig. S2b-e online). Overall, transcriptome data of 4795 protein-coding genes from 82 cells was retained (Fig. S2f  online). IL2RA and ACTB exhibited unimodal expression, with FOXP3 and CTLA4 having bimodal distribution, suggesting phenotypic variation (Fig. 1a). IL2RA expression was unimodal among Treg subpopulations, which suggests CD25 can not be used as a surface marker to distinguish stable Treg and unstable Treg.
We performed dimensional reduction on the gene expression matrix using the Rtsne package, showing the formation of two clusters (Fig. 2b). To explore the inherent difference among these Treg subpopulations, we carried out differential expression analysis and hierarchical clustering, finding two stable clusters each with their unique signature genes (Fig. 2c). Unbiased trajectory ordered cells according to Treg functionality (Fig. 2d). According to FOXP3 and Treg functional gene expression levels, we defined these Treg subpopulations as stable Treg, unstable 1 and unstable 2 (Fig. 2e). Unstable Treg cells decreased Treg signature gene expression including FOXP3, IL2RA, CTLA4, TNFRSF18, and Helios (Fig. 2e). The unstable 1 cluster was transitional, which shared partial expression profiles with both the stable and unstable 2 clusters (Fig. 2c). t-SNE projection also showed an increase in the correlation between the unstable 1 and unstable 2 populations suggesting that unstable Treg cells are a more homogeneous population compared with stable Treg cells.

TIGIT + Treg cells are more stable under IL-6 stimulation
It was important to identify an appropriate cell surface biomarker to distinguish between stable and unstable Treg cells. Based on single-cell sequencing data, we found a series of genes that correlated with FOXP3 expression, of which TIGIT was found expressed at a low level in unstable Treg cells (Figs. 2e, 3a). The profile of human cord blood naïve Treg cells showed that TIGIT + Treg cells compared to their TIGIT À counterparts had no significant difference in the expression of Treg signature genes (Fig. 3b and c). However, TIGIT À Treg cells significantly downregulated Treg signature gene expression after anti-CD3/CD28, IL-2 plus IL-6 treatment ( Fig. 3d and e). In order to investigate the source of TIGIT À Treg, we sorted primary TIGIT + and TIGIT À Treg cells, then treated with anti-CD3/CD28, IL-2 plus IL-6. We found that the source of TIGIT À Tregs subpopulation mainly derived from primary TIGIT À Treg cells. TIGIT À Treg cells significantly downregulated Treg signature gene expression (Fig. 3f). FOXP3 as a crucial transcription factor can bind many Treg signature genes locus and influence their expression including TIGIT, CTLA4 and IL2RA [39]. The absent of FOXP3 directly decreased the transcription of Treg signature genes in TIGIT À Treg cells. These results indicate that TIGIT + Treg cells are more stable and have more powerful suppressive function.

Unstable Treg cells have increased cell proliferation and proinflammatory cytokine expression
Unbiased gene-set enrichment analysis showed that cell proliferation-associated gene-sets were up-regulated in unstable Treg cells (Fig. 4a). Furthermore, unstable Treg cells increased Ki-67 expression compared with their stable counterparts (Fig. 4b). Increased glycolysis can provide biosynthetic intermediates for unstable Treg cells growth and proliferation (Fig. 4a) [40]. In addition, unstable Treg cells decreased expression of Treg immunosuppressive factors (IL-10, GZMB and PRF1) and upregulated pro-inflammatory cytokines (IL-15, IL-16 and IL-18) ( Fig. S3a, b online). Unstable Treg cells also highly expressed TNFa and IL-17A (Fig. 4c, d). These results indicated that unstable Treg cells lose Treg signature and gain a Th17-like phenotype after treatment with IL-6.

Stable Treg cells have elevated expression of receptors and signaling toward inflammation signals
Next, we sought to investigate the transcriptional profile of stable Treg populations. Stable Treg cells highly expressed coactivation markers including CD48, CD44, CD28 and PD1, indicating that the stable cells were activated (Fig. S3c online). Furthermore, we found stable Treg populations also upregulated CCR4 expression. CCR4 high expressed Treg cell represented a highly immunosuppressive subset (Fig. S3d online) [41]. Gene-set enrichment analysis showed that the most significant functional gene classes with differences were cytokine-associated pathways between stable and unstable Treg cells (Fig. 5a). The IL-2/IL-2R signal has been known to a crucial for Treg homeostasis in the periphery. Stable Treg cells had upregulated IL-2-STAT5 signaling that could maintain Treg stability. Moreover, stable Treg cells displayed higher levels of inflammatory cytokine signaling pathways compared to unstable Treg cells, including IFN-c, IL-6 and TNF-a (Fig. 5a). We assessed expression of IL6R, IL1R1, INFGR1 and TNFRSF1A by flow cytometry and found that expressions of all of these receptors were increased in stable Treg (Fig. 5b).

Co-expression network analysis reveals CYP1A1 as a critical regulator of Treg cell function
Single cell RNA-seq provides an approach to deep understand any potential candidates for modulating human Treg cell function. We selected high expression genes (frequency > 20%, FPKM > 5) to investigate novel regulators of Treg cells. To unveil the coexpression relationships of genes, we performed weighted gene co-expression network analysis using the genes filtered from above [42]. In total, 4795 genes were selected for co-expression analysis. We found a key transcript module that co-varies with known regulatory genes including FOXP3, IL2RA, CTLA4, BATF and IL10 (Fig. 6a). We prioritized the candidates with a computational ranking scheme. Based on ranking, we chose CYP1A1 for further investigation. CYP1A1 have a higher expression level in Treg cells compared with conventional T cells (Fig. 6b). The strong correlation between the expression of CYP1A1 with FOXP3 and TIGIT was confirmed at the protein level (Fig. 6c). IL-6-treated Treg cells showed that nearly all stable Treg cells expressed CYP1A1. These data indicate that CYP1A1 may be a regulator of Treg cell stability and function.

CYP1A1 is an important regulator of human Treg cell function
Cytochrome P450 (CYP) proteins are heme-thiolate enzymes involved in multiple biosynthesis and metabolic processes. CYP1A1 is a member of the CYP superfamily. To assess the role of CYP1A1 in Treg cell function, we knocked-down CYP1A1 with shRNA-bearing lentivirus (Fig. 7a). The expression of regulatory genes including FOXP3, CTLA4 and IL2RA were substantially downregulated both at the RNA (Fig. 7b) and protein levels (Fig. 7c, d). To examine the immunosuppressive capability of CYP1A1-depleted Treg cells, we carried out a suppression assay. Results indicated that CYP1A1-depleted Treg cells displayed impaired immunosuppressive activity to effector T cell (Fig. 7e). Furthermore, we investigated the function of CYP1A1 in Treg cells using the CYP1A1specific inhibitor Alizarin. We assessed expression of FOXP3, CD25 and CTLA4 by flow cytometry and found that these genes were downregulated in a dose dependent manner, indicating that CYP1A1 plays a crucial role in Treg cells (Fig. 7f). Previous findings have shown that unstable Treg cells can lose FOXP3 expression and produce higher amounts of inflammatory cytokines. We questioned whether the knockdown of CYP1A1 could promote the expression of inflammatory cytokines in Treg cells. We observed that IL17A expression increased in CYP1A1-depleted Treg cells with IL-6 treatment (Fig. 7g). Furthermore, there are no IL17A production in CYP1A1-depleted Treg cells without IL-6 treatment which indicated CYP1A1 can't divert Treg cells to Th17 cells without IL-6 stimulation (Fig. 7g). These results collectively suggest that CYP1A1 maintains Treg cell suppressive function and contributes to Treg cell trans-differentiation into Th17-like cells in an IL-6-dependent manner.

Discussion and conclusion
Single-cell RNA-seq has been successfully utilized to study cellular heterogeneity in several tissues at the genomic and/or transcriptomic level. In this study, we utilized single-cell RNA-seq to investigate the heterogeneity of Treg after IL-6 treatment. Covariation analysis highlighted novel regulators of Treg cell function and stability, which have not been detected by previous population level approaches [30].
Umbilical cord blood immune cells are mainly immature and less active against antigen responses partly due to the abundance of regulatory T cells [43]. We observed that Treg cells from umbilical cord blood have enhanced immunosuppressive capability and proliferation compared with adult peripheral Treg cells. Therefore, we took advantage of naïve Treg cells derived from human umbilical cord blood to dissect transcriptome variation of IL-6stimulated Treg cells at the single cell level.
Single-cell RNA sequencing provides transcriptome data in individual Treg cells. Our data demonstrates that in the presence of high concentrations of IL-6, ''unstable" Treg cells lost their suppressive capacity and expanded in number. However, ''stable" Treg cells retained suppressive activity and had higher expression of responsive elements to inflammatory cytokines. We speculate that unstable Treg cells may contribute to fighting pathogenic disease, while stable Tregs have higher sensitivity to and respond to inflammatory cues to migrate to inflammatory sites to resolve inflammation. These results reveal an understanding of Tregs and their microenvironment and the potential dynamics of Treg homeostasis during inflammation. These findings also reveal an underlying mechanism for future immunotherapies aimed at fine-tuning Treg-mediated immune responses, including cancer and autoimmune diseases. TIGIT defines a highly suppressive and stable Treg cell population. TIGIT + Treg cells are more stable under inflammatory conditions. We observed that the majority of IL-6 stimulated TIGIT À Tregs were derived from TIGIT À Treg cells and displaying a Th17 phenotype, while only a small proportion of TIGIT À Tregs could also be derived from TIGIT + stable Treg cells. We found that the specialized TIGIT + and TIGIT À Treg cell subsets display a certain degree of plasticity, and FOXP3 as a crucial transcription factor to maintain Treg cell stability is involved in this progress. IL-6 can balance RORct and FOXP3 expression in TIGIT À Tregs by promoting FOXP3 protein degradation and relieving FOXP3 + Treg mediated immune suppression. In addition, TIGIT + Tregs were found to express very high levels of immunosuppress cell surface molecules including TIGIT, CTLA4 and IL2RA, which are key mediators of Treg suppression in vivo. Our results indicate that TIGIT + Treg cells may serve as desirable Treg cell platform for further synthetic biology based cell modifications to make CAR-Tregs for treating graft-versus host disease (GVHD) as well as other major autoimmune diseases [44,45].
Single-cell RNA sequencing can identify novel gene regulators which are hard to be found in traditional bulk RNA sequencing. The transcription factor FOXP3 can drive Treg lineage genes expression and determine the immunosuppressive capability of Treg cells [7,46]. In our model, CYP1A1 is strongly co-expressed with known Treg regulatory genes. Consistently, the lack of CYP1A1 significantly impaired Treg cell function. Interestingly, it has been recently found that the impaired expression of CYP1A1 in mice can exhaust aryl hydrocarbon receptor (AHR) ligands [47]. AHR is a liganddependent transcription factor. Many substrates for CYP1 enzymes are AHR ligands. AHR is located in the cytosol when the cells are lack of its ligands. Ligand binding can activate AHR to initiate its transcription activity on multiple AHR-responsive genes including CYP1A1. Persistent CYP1A1 enzyme reaction is crucial for AHR activation as a positive feedback loop for the regulation of AHR target genes. Meanwhile, AHR can control Treg cell generation and FOXP3 expression, indicating that CYP1A1 may regulate Treg cell function and proliferation in an AHR-dependent manner. CYP1A1 catalyze many reactions involved in eicosanoid, cholesterol, sterol, lipid and steroid synthesis and metabolism. The expression of CYP11A1 is increased in Th2 cells, and published data has also shown that Th2 cells are capable of de novo steroid synthesis by CYP11A1 [29]. Therefore, CYP1A1 may be involved in steroid biosynthesis and contributes to Treg function and proliferation. Our data suggested that CYP1A1 and TIGIT could be considered as a potential target for immunotherapy of autoimmune diseases and cancer. Future studies should test the functional role of CYP1A1 in regulating Treg cell stability and functions in vivo by using Treg specific CYP1A1 conditional knockout mice models under disease settings.
In summary, we used single-cell RNA sequencing to investigate the heterogeneity and plasticity of Treg cells during inflammation using IL-6 as a stimulus. We defined and functionally analyzed two distinct Treg subpopulations, TIGIT + ''stable" Treg and TIGIT À ''unstable" Treg cells. Meanwhile, we found CYP1A1 is novel regulator of Treg cell function by single cell transcriptome analysis. Therefore, targeting CYP1A1 might have broad potential and application in diseases that require the perturbation of Treg function.  Research in Dr. Li's laboratory mainly focuses on understanding molecular mechanism underlying the functional stability, plasticity and balance of FOXP3 + Regulatory T cells (Treg).