Single-Molecule DNA Methylation Reveals Unique Epigenetic Identity Profiles of T Helper Cells

In general, adaptive immunity in both protective and autoimmune settings is orchestrated by the emergence of polarized effector CD4 T helper (Th) lymphocyte subsets of specialized function. Naive CD4 Th0 cells under the influence of cytokines and TCR signaling can give rise to different effector cells including Th1, Th2, Th17, and T regulatory (Treg) cells (1, 2). Furthermore, plasticity is an important property of Th cells, whereby they adapt their functions in response to a changing environment forming a continuum of polarized phenotypes (3). Hence, cells with intermediate Th1/Th17 and Th17/Treg phenotypes have been described in both mice and humans in association with autoimmune diseases and protective functions, respectively (4).

Both identity and plasticity of Th cells are regulated in part by epigenetic mechanisms (2, 5). Epigenetics refers to the heritable layer of information on top of a DNA sequence that regulates gene expression and makes each cell type unique. Therefore, mapping epigenetic modifications represents an ideal way to discriminate cell types by their “epigenotype” (6, 7). DNA methylation (5mC) is the most stable and most widely studied epigenetic mark and is involved in gene silencing (8). DNA hydroxymethylation (5hmC) is relatively understudied but recently has been linked to gene activity during differentiation (9, 10). Because of their relative stability (11) and cell-type specificity (9), mapping 5mC/5hmC landscapes has the power to improve our understanding of Th cell identity and plasticity and aid in the development of novel biomarkers for these cells.

Despite the evidence for 5mC/5hmC as a marker of cellular identity, in particular in Th cells, technical limitations have prevented this knowledge from reaching clinical use. Simultaneous detection of these cytosine modifications in CpG context (5mCpG and 5hmCpG) has been made recently available by the next generation of nanopore long-read sequencers, which are able to provide such information in native DNA. In this article, we take advantage of nanopore sequencing to simultaneously profile 5mCpG and 5hmCpG in native DNA on naive Th cells polarized under Th0, Th1, Th2, Th17, and Treg conditions. In addition, given that transforming growth factor β (TGF-β) is known to prevent pathogenic Th17 cell differentiation (12), we performed targeted sequencing on both in vitro– and in vivo–generated Th17 cells escaping TGF-β control. This study opens new avenues toward the implementation of 5mC/5hmC in clinical immune profiling and identification of pathogenic immune cell subsets.

Foxp3-GFP reporter male mice (13) between the ages of 6 and 12 wk of age were used for in vitro Th polarization. IL-17a–CRE mice (14) were bred with mice bearing floxed alleles of TgfbR2 (15). These so-called TGFbR-KO mice and their TGFbR-WT counterparts (littermates without floxed alleles) were then crossed with Rosa-26-stopFlx/Flx‐yfp mice to enable tracking by monitoring yellow fluorescence protein (YFP) expression, which represents the cells that either are expressing or expressed IL-17a.

All animals were housed in ventilated racks in pathogen-free facilities at the Institut de Génétique Moléculaire de Montpellier (Montpellier, France) and Plateforme Petit Animal du Cancer Research Center of Lyon (Lyon, France). The animal care and experiments were approved by the local and national animal facility institutional review boards in accordance with French national and ARRIVE guidelines (CLB 2017-017).

Small intestine colon and mesenteric lymph nodes from TGFbR-WT and TGFbR-KO mice were dissected. Fat was removed, and the intestines were longitudinally opened and washed in 1× PBS (Life Technologies). The intestines were cut into small pieces and incubated with 5 mM EDTA, 1 mM DTT (D0632, Sigma–Aldrich) at 37°C. Epithelial cells were then separated from intraepithelial lymphocytes with a 44%/67% Percoll gradient (P1644, Sigma–Aldrich) run for 20 min at 1,300 × g. The tissues were then digested in RPMI medium (Life Technologies) containing 20% FBS (Life Technologies), DNAse I (Roche) at 100 µg/ml and collagenase from Clostridium histolyticum (C2674, Sigma–Aldrich) at 1 mg/ml. Intestinal lamina propria lymphocytes were then separated on a 44%/67% Percoll gradient run for 20 min at 1,300 × g. Mesenteric lymph node cells were ground between wire mesh. CD4⁺ T cells were enriched using the CD4⁺ T cell isolation kit (Miltenyi Biotec) and then labeled with CD4 PE (GK1.5, eBioscience) TCR-β APC (H57-597, BD Biosciences), CD45 APC-Cy7 (30-F11, BD Biosciences), and DAPI (Eurobio Scientific). The cells were purified on ARIA II (BD Biosciences). For in vitro polarization, CD4⁺ T cells were purified from Foxp3-GFP mice using the MACS CD4⁺ T cell negative selection kit (Miltenyi Biotec), and I CD4⁺ T cells from mice were then sorted on the basis of a CD4⁺CD8⁻CD62L⁺CD44⁻GFP⁻CD25⁻ expression profile on a FACSAria flow cytometer (BD Biosciences).

Purified naive CD4⁺ T cells from male Foxp3-GFP mice were activated under different polarizing conditions for 4 d: Th0 (neutral), Th1, Th2, Th17, and Treg, as well as Th17s differentiated without TGF-β1 (R&D Systems, catalog no. 240-B) (Th17–no TGF-β) representing an intermediate cell phenotype (Supplemental Tables I, II). Activation was performed using plate-bound anti-CD3 (clone 145-2C11, 1 μg/ml; BioXcell, catalog no. BE0001-1) and anti-CD28 (clone PV-1, 1 μg/ml; BioXcell, catalog no. BE0015-5) monoclonal Abs in RPMI 1640 medium (Life Technologies) supplemented with 10% FCS, 1% penicillin/streptomycin (Life Technologies) and β-mercaptoethanol (50 μM). The cells were split 3 d later with medium supplemented with rhIL-2 (100 U/ml). The cells were maintained in a standard tissue culture incubator containing atmospheric O₂ and 5% CO₂. 5mCpG raw data from Tregs included in this study were previously published by us under the exact same experimental conditions (5). These data were reanalyzed to compare an additional CD4 Th subset using updated basecalling algorithms and analysis pipelines.

Immunophenotyping of cells was performed with fluorochrome-conjugated Abs, and intracellular staining was performed after fixation and permeabilization (ThermoFisher). The cells were labeled with fixable viability dye prior to fixation, followed by at least 3 h of intracellular staining with the following anti-transcription factor Abs: Foxp3-PeCy7, T-Bet-PE, Rorɣt-BV650, and Gata3-APC (Supplemental Table I). Cytokine production (IFN-γ, IL-4, IL-17A, and IL-10) was also assessed by a cytometric bead array kit (BD Biosciences) on supernatants collected at day 3 of polarization. Data analysis was performed using FlowJo (Tree Star software) and FCAPArray software (cytometric bead array analysis).

DNA was extracted from 1 × 10⁶ CD4 Th cells/replicate using the Qiagen DNA mini kit according to the manufacturer’s instructions (Qiagen, Hilden, Germany). DNA was extracted from 10 ng of mouse brain tissue as previously described (16) using NEB Monarch high m.w. DNA extraction kit according to the manufacturer’s instructions (NEB catalog no. T3050). Viable CD45⁺CD4⁺YFP⁺ cells were sorted from the small intestine, lamina propria, Peyer’s patches, and mesenteric lymph nodes of 6- to 9-mo-old mice using a FACS ARIA II (BD Biosciences). A Beckman GenFind V3 DNA extraction kit (Beckman Coulter) was used to isolate DNA from 120,000 TGFbR-WT and 180,000 TGFbR-KO sorted cells, following the manufacturer’s instructions.

DNA was sheared using g-tubes (Covaris, USA) to ∼10 kb before library preparation with ligation sequencing kit (ONT catalog no. SQK-LSK109) followed by whole-genome sequencing on nanopore MinION flow cells with pore 9.4.1 chemistry for 72 h. Raw read quality was determined with pycoQC (17).

A double cutting approach was used whereby 2× Streptococcus pyogenes Cas9 Alt-R trans-activating CRISPR RNAs (crRNAs) were designed for both upstream and downstream of each target loci as previously described (5, 18) (crRNA details and sequences are provided in Supplemental Table III). crRNAs were pooled in equimolar ratios (100 µM) and annealed to S. pyogenes Cas9 Alt-R tracrRNA (100 µM) (Integrated DNA Technologies, Coralville, Iowa). crRNA·tracrRNA pool (10 µM) was incubated with Cut smart buffer (New England Biosciences catalog no. B7204) and Alt-R S. pyogenes HiFi Cas9 nuclease V3 (62 µM) forming crRNA-tracrRNA-Cas9 ribonucleoprotein (RNP) complexes. High m.w. DNA (3,000 ng) was prepared by blocking available DNA ends with calf intestinal phosphatase (NEB catalog no. M0525). DNA and dA tails on all available DNA ends were cleaved by RNPs and Taq polymerase (New England Biosciences catalog no. M0273) and dATP (New England Biosciences catalog no. N0440), activating the Cas9 cut sites for ligation. Nanopore native barcodes (Oxford Nanopore Technology catalog no. EXP-NBD104) were ligated to available DNA ends with NEB blunt/TA ligase master mix (New England Biosciences catalog no. M0367), and up to four samples were pooled before ligation of nanopore adapters with NEBNext Quick T4 DNA Ligase (New England Biosciences catalog no. E6057) and loading onto a primed MinION flow cells (pore chemistry 9.4.1) and sequencing for 72 h. Raw read quality was determined with pycoQC (17).

For adaptive sequencing we used the “UNCALLED” methodology in targeted reference genes (19). To this end, a reference gene panel was developed based on known targets of TGF-β (20). Gene sequences were retrieved from the mouse mm10 assembly, including a 10-kb flanking region up and downstream of each region of interest, and masking repeated elements. This method did not require Cas9-based library preparation. Instead, during sequencing, UNCALLED is able to rapidly match streaming of nanopore current signals to the reference sequence, enriching in this way for our selected 139 targets. A full reference list and the UNCALLED loci are listed in Supplemental Table IV.

Raw fast5 files were converted into pod5 format and basecalled with Dorado (version 0.4.2, ONT), allowing the simultaneous delineation of 5mCpG and 5hmCpG using a high-accuracy model ([email protected], ONT) for the detection of modified bases (18). Dorado was also used for simultaneous alignment to the mm10 assembly and demultiplexing, followed by sorting and indexing with Samtools (version 1.13). Modkit (version 0.2.2) was used for piling of methylation data into bed files.

Data processing and statistical analyses were performed using R/Bioconductor (R version 4.3.1). Bed files containing modified cytosine information were transformed into a BSseq object for differential methylation analysis with dispersion shrinkage for sequencing data as previously described (5, 21). Tests for dispersion shrinkage for sequencing for differential methylation at single CpG sites were evaluated using a Wald test on the coefficients of a β-binomial regression of count data with an “arcsine” link function. Differentially methylated regions (DMRs) were defined as those loci with at least three CpG sites within a distance of less than 50 bp and with changes in >50% of all CpG sites exhibiting p value < 0.05. DMRs were plotted using the plotDMRs function of the dmrseq R package (22, 23). Off-targets were analyzed for transcription factor (TF) activity using the cistrome information associated with the MIRA R package (24). After extracting all genomic binding regions corresponding to selected TFs, MIRA functions were used for aggregation of DNA methylation data, estimation of TF scores, and visualization of inferred TF activity. Only TFs described as high confidence were used, with a coverage of at least 100× and at least 10,000 annotated targets.

To study the relationship between DNA methylation and classic Th cell identity, we polarized naive cells toward various Th fates in vitro. Importantly, this approach does not lead to a terminal differentiation state but rather a polarization toward different fates (25). Naive cells were isolated by FACS sorting and activated under polarizing conditions toward the Th0, Th1, Th2, Th17, and Treg phenotypes using cytokine cocktails. To test the sensitivity of this approach to detect potential pathogenic Th cell phenotypes, naive T cells were also polarized under Th17 conditions in the absence of TGF-β, referred to as “Th17–no TGF-β” cells (Fig. 1A).

Specific transcription factor staining (Fig. 1B, 1C) and cytokine production (Supplemental Fig. 1A) associated to Th cells were performed to estimate the efficacy of polarization for each condition. Expression patterns of transcription factors T-bet, Gata3, and Foxp3 appeared selectively upregulated in Th1, Th2, and Treg polarizing cells, respectively (Fig. 1B, 1C), and correlated with the cytokine signature detected (IFN-γ and IL-4 for Th1 and Th2 conditions, respectively; Supplemental Fig. 1). Expression of Rorgt was consistently upregulated in Th17 and Th17–no TGF-β polarized cells (Fig. 1B, 1C). Mean fluorescence intensity of Rorgt increased by 1.8-fold in Th17 and 1.7-fold in Th17-no TGF-β as compared with control conditions. Furthermore, only cells polarized toward a Th17 fate produce significant amounts of IL-17 and did not produce IFN-γ or IL-4 at a detectable level (Fig. 1C). Taken together, these data indicate an efficient polarization of Th17 and Th17–no TGF-β conditions.

To characterize the DNA methylation landscape of polarized Th cells, we developed a Cas9-guided capture method coupled with nanopore sequencing (Fig. 1A). We manually selected key genes that could potentially be used to classify different Th subsets. This panel included DNA regions in which differential epigenetic modifications in Th cell subsets have been previously reported (26–28). In general, these genes code for major TFs and cytokines specific to Th cell subsets, namely IL-17a, Maf, Rorgt/Rorc, and Zfp362 (28) (Th17); Foxp3, IL-10, and Itgb8 (Treg); Gata3, IL-4, and Il10 (Th2); and Tbx21/T-bet and IFN-γ (Th1). Cas9-guided enrichment captures gene length regions prior to single-molecule native DNA sequencing with nanopore sequencing. This involved the design of two sgRNAs upstream and downstream of each gene of interest employing a double cutting approach. It was then followed by synthesis of Cas9-guided RNPs to enrich for targets prior to nanopore sequencing (Fig. 1D) as previously described (5, 7, 21). Up to four samples from each polarization assay (i.e., Th0, Th1, Th2, Th17, Th17–no TGF-β, and Tregs) were barcoded and multiplexed on each MinION flow cell, and profiled for 5mCpG and 5hmCpG landscape of whole gene length single DNA molecules. In line with previous work on other cell types, the coverage of 10× was sufficient to detect a 25% change in DNA methylation (7) and was therefore set as the minimum target coverage for this study. The typical sequencing yield achieved for each of these nanopore runs was ∼1 GB, which is in accordance with the literature (7, 18). Importantly, we exceeded the minimum coverage threshold of 10× for target regions (average coverage, 45×; Fig. 1D, Supplemental Fig. 1B) and obtained full gene length reads (N₅₀ = 21 kb) on all targets (Fig. 1E, Supplemental Fig. 1C); this was sufficient for further processing to detect differentially methylated regions. These data highlight the efficacy of Cas9 targeted sequencing for capturing important regions for Th cells.

To basecall raw reads and detect modified bases, we used the latest high-performing nanopore tool, Dorado (https://github.com/nanoporetech/dorado), combined with high-accuracy models to call DNA modifications. Dorado works by extracting modified base calls from raw nanopore reads and anchoring the information-rich basecalling neural network output to a reference genome (mm10). Global patterns of methylation were compared between different nanopore sequencing runs, and no differences were detected (Fig. 1F, 1G). In all sequencing runs, 5mCpG displayed a typical bimodal distribution with peaks close to zero and 100% methylation and a variable fraction of CpG sites with intermediate levels of methylation (Fig. 1F). In contrast, 5hmCpG was characterized by very low levels in all sequencing runs, with a much lower fraction of detectable methylated sites above 10% (Fig. 1G).

Differentially methylated loci (DMLs) between subsets were defined as individual CpG sites with significant change (p value < 0.05) of at least 10% in any direction. The number of distinct DMLs identified in Th1, Th2, Th17, and Treg cells when comparing to Th0 cells, were 305, 359, 428, and 490, respectively. Furthermore, a comparison of Th17 cells polarized in the presence or absence of TGF-β (i.e., Th17–no TGF-β cells) identified 378 DMLs, highlighting a potential effect of TGF-β alone on modulating DNA methylation during T cell polarization.

To obtain a more robust measure of 5mCpG variation, we searched for DMRs, defined as genomic ranges of at least 50 bp and a minimum of 3 CpG sites, where at least half of the CpG sites are significantly different (p < 0.05 and 10% methylation change) between two cell subsets. This analysis resulted in 20 DMRs that distinguished Th0 cells from the different polarized subsets (Supplemental Table V). Two additional DMRs were identified in Th17–no TGF-β cells relative to Th0 cells and Th17 cells polarized in the presence of TGF-β. The length of these DMRs averaged 214 bp and ranged from an off-target 53-bp region in the Tnk2 locus to a large 906-bp DMR in the Il4 locus displaying 60% demethylation in Th2 cells relative to Th0 cells. In terms of CpG content, they ranged from 4 to 15 sites and an average of 6 CpG sites per DMR. In line with the selection of the targeted DNA regions, DMRs were largely negative compared with Th0 cells, highlighting that in these identity genes, 5mCpG decreased during T cell differentiation in correlation with increased gene expression. Significantly, hierarchical clustering of the resulting 22 DMRs was sufficient to cluster Th cells into their respective subsets (Fig. 2A). Interestingly, Th2 cells were the most divergent of the different cell subsets in this clustering analysis (Fig. 2A).

Plots of each DMR (Fig. 2B–F) illustrate the specificity of methylation profiles for each cell type. Indeed, DMRs upstream of the Foxp3 promoter and in a CpG island (CGI) shore on Tbx21 were unique to Treg and Th1 cells, respectively (Fig. 2B, 2C). Eight DMRs distinguished Th2 cells, three of them in the Il4 gene body and its promoter (Fig. 2D, upper panel) and the others in Il10, Gata3, and close to the Foxp3 locus. Six DMRs were associated with Th17 polarization, including one in a distal enhancer close to Zfp362 and two in the body of the Rorc gene (Fig. 2E). The same Zfp362 enhancer distinguished Th17 cells polarized in the presence or absence of TGF-β (Fig. 2F), with Zfp362 remaining hypomethylated in Th17–no TGF-β cells at levels close to Th0 cells.

These data highlight the power of characterizing DNA methylation profiles in a small subset of genes as a means to determining Th cell identity. They also suggest that 5mCpG profiles can be used to evaluate the specific impact of a cytokine on Th differentiation.

Although 5hmCpG is generally associated with increased gene expression when present in coding regions (21), the relationship between 5mCpG and 5hmCpG in CD4 Th gene regulation is poorly described. This is of importance because 5mC and 5hmC have been shown as key players in tissue differentiation and cellular plasticity (9). Therefore, we took advantage of nanopore sequencing to simultaneously detect 5mCpG and 5hmCpG on the same native nanopore reads, as described in Materials and Methods. Although the detection of 5mCpG with nanopore sequencing has been validated, the detection of 5hmCpG has not. Given that tissue derived from the CNS is known for being high in DNA hydroxymethylation (29, 30), we used data obtained from DNA extracted from the brain of P13 OF1 mice as a positive control for 5hmCpG. Using the same pipeline described above for calling modified bases, a small fraction of CpG sites (3%) displayed 5hmCpG of at least 10% in brain tissue. This percentage increased to 4.2% when aggregating the 5hmCpG data corresponding to forebrain enhancer regions. This was significantly higher than the 3.3% average 5hmCpG identified at control adipose enhancer regions (p < 0.05; Supplemental Fig. 1D, 1E). Although this is far from a rigorous validation of nanopore sequencing–based 5hmCpG accuracy, we proceeded to determine the 5hmCpG landscape of CD4 Th cell subsets in an exploratory manner.

No DMRs (5hmCpG) were identified with the same criteria used earlier for 5mCpG. Despite this, 42 regions contained significant CpG sites in discrete regions, regardless of their magnitude of methylation change (Supplemental Table VI). Together, these regions do not discriminate cell subtypes, as opposed to 5mCpG (Fig. 3A, 3B). Again, two of the three Th2-polarized samples are clustered apart, in particular due to few loci mapping to Il4, Il10, Zfp263, and Gata3 (Fig. 3A). Of note, most of these regions displayed higher 5hmCpG levels in Th2 cells, in contrast to their lower 5mCpG levels described above. Based on these results, 5hmCpG does not add discriminatory capacity to 5mCpG at the loci and conditions of this analysis.

In addition to targeted regions, nanopore sequencing captures ∼1 GB of off-target reads in each sequencing run. To strengthen the identification of Th cell subsets using DNA methylation, we determined the DNA methylation status of TF-binding sites in on- and off-target reads with Bioconductor package MIRA (methylation-based inference of regulation activity) (24). MIRA aggregates the DNA methylation values in all known binding sites for chosen TFs and uses this profile to produce a score. The score infers the level of regulatory activity on the basis of the shape of the DNA methylation profile, based on the observation that DNA methylation tends to be lower in regions where TFs are bound (31). In a first step, TFs were chosen matching the known transcriptome orchestrators of Th cell subtypes Tbx21 and Stat1 (Th1), Gata3 (Th2), Stat3 and Runx3 (Th17), and Foxp3 (Treg) (32, 33). Bcl6 is a TF for T follicular helper cells (33) and was used as a negative control for target Th cells. We aggregated 5mCpG located at known binding sites for the chosen TFs, representing thousands of regions (Fig. 4, left panels). In general, we observed a hypomethylation dip for all TFs, indicating regulatory activity. In addition to plotting MIRA profiles, we determined the MIRA scores for each cell subset presented as boxplots next to the aggregated methylation plots (Fig. 4). Three of the selected TFs fell below statistical significance (i.e., FOXP3, GATA3, and RUNX3), with Tregs displaying the highest MIRA scores for FOXP3 binding sites; Th2, Tregs, and Th17 displaying high scores for GATA3 binding sites; and Th17 cells having the highest score for RUNX3.

To provide an unbiased analysis of off-target sequences, we next performed the same analysis for all high-confidence TFs available in the MIRA database. Of 877 TFs, only 181 are classified as high confidence by MIRA’s own classification. After keeping only those TFs with more than 10,000 targets across the genome and at least 100× coverage in our data, we inferred the scores for 87 TFs. Only two TFs were significantly different after p value correction (false discovery rate < 0.05): IKZF1 and PAX5 (Fig. 4). In both cases, Th17 cells displayed the highest MIRA scores. Together, these data show that the methylation status of TF binding sites can be used as an indicator of the activity of Th cell subsets, further strengthening the relationship between 5mC and Th cell identity.

Our results above suggest that TGF-β is partially involved in shaping the methylome landscape during polarization of Th17 cells. In order to validate our results in vivo, we used Th17 cells present in the gut of Il17aCRE;Tgfbr2 fl/fl Rosa-26 stop fl/fl⁻ yfp (TGFbR-KO) mice and Il17aCRE;Rosa-26 stop fl/fl⁻ yfp (TGFbR-WT) mice. These cells escape TGF-β signaling control and can be purified from the tissues based on YFP expression. Their methylation profiles were analyzed through a modified method of targeted sequencing in which single molecules are selected during nanopore sequencing with coupled real-time basecalling (see details in Materials and Methods). This so-called adaptive sampling was used to extend the selection of genomic regions from the original panel of 10 genes to 139 annotated genes, selected because of their known role in TGF-β signaling (20). We compared 5mCpG in TGFbR-KO Th17 cells to their WT counterpart. A total of 330 DMRs distinguished the two conditions, with 21 of them within the original 10-gene panel (Fig. 5A; Supplemental Table VII). In line with the pathogenic potential linked to the absence of TGF-β control on Th17 cells (12), Tbx21, encoding for T-bet, was the top most significant of these 10 genes, with 8 DMRs falling within its locus and 7 of them showing hypomethylation in TGFbR-KO cells (Fig. 5B). Moreover, two of these DMRs displayed differences in methylation of more than 60% (hypomethylated in TGFbR-KO cells relative to TGFbR-WT cells). Similarly, in line with a epigenetic signature of in vivo generated Th17 cells (28), three DMRs were identified in the Zfp362 locus, including a large hypomethylated region upstream of the Zfp362 promoter that overlaps with the DMR identified in in vitro differentiated Th17 cells in the absence of TGF-β (Fig. 5B). One of these DMRs matched a distal enhancer region and was at least 60% less methylated in TGFbR-KO cells relative to TGFbR-WT cells.

CCCTC-binding factor5hmCpG was also studied in the Th17 cells from the mouse gut. Overall, 5hmCpG was detectable in a higher proportion of CpG sites and at higher levels than in vitro polarized cells, although no DMRs were detected in any of the genes of the original 10-gene panel. When expanding the search to the 139 targeted regions, three 5hmCpG DMRs were identified in two different genes: Zeb2 and Pparg. One of these, a 173-bp region in the Zeb2 locus, spanned 8 CpG sites with 5hmCpG levels that increased from an average of 3% in TGFbR-WT cells to 38% in TGFbR-KO cells. The same region was identified as a 5mCpG DMR hypomethylated in TGFbR-KO cells (Supplemental Table VII). Thus, 5hmCpG can indeed be detected in Th17 cells at discrete loci in vivo. This set of data highlights the potential to exploit native DNA methylation profiling to study physiological and pathological features of Th17 cells present in tissues.

DNA cytosine modifications are known to be cell type specific (1). However, the role of these modifications in Th cell identity is not well understood. The current study addressed the relationship between DNA methylation and Th cell identity, characterizing 5mCpG and 5hmCpG of murine CD4 Th cell subsets (Th0, Th1, Th2, Th17, Th17–no TGF-β, and Tregs) by long read single-molecule native DNA sequencing. We propose a novel approach able to identify differentially methylated regions expected for each Th cell subset. Using this method, our data revealed that 5mCpG profiles from different loci provide a map specific to each in vitro–polarized Th cell subtype. Moreover, we confirmed the efficiency of nanopore targeted sequencing to discriminate pathogenic and nonpathogenic Th17 cells in vivo.

We developed a Cas9 targeted sequencing to enrich important loci for Th cell identity. This method obtained high coverage of target genes compared with short read sequencing approaches, as well as whole gene length reads. Cas9 enrichment is an innovative technique that has been adopted to answer a variety of biological questions ranging from studying epigenetic heterogeneity in viruses (7), important loci for cancer (18), and neuromuscular disorders (34). Regardless of the question, Cas9 targeted enrichment coupled with nanopore sequencing obtains anywhere from 0.5 to 2 GB of data and tenths to thousandths of on-target reads. In the current study, we observed variability between coverage of different target loci achieving the highest coverage of Il4 and the lowest of Itgb8. This coverage disparity is likely due to differences in size of target regions. Il4 is a smaller gene, and thus to obtain all regulatory information nearby we aimed to enrich in ∼7 kb, whereas Itgb8 is a larger gene and thus required enrichment of >100 kb. Shorter fragments elicit a higher coverage due to the nature of nanopore sequencing’s preferentially processing shorter reads. Nevertheless, the on-target reads obtained in the current study fall well within that of the literature (18, 34).

Our data show that 5mCpG in a selection of genes can accurately classify polarized Th cell subsets. Importantly, capturing long reads with high coverage of important loci revealed that DMRs were not only specific for a single subset; instead, some regions exhibited a range of methylation profiles. Other techniques relying on sodium bisulfite conversion are not able to distinguish 5mC from other modified bases such as 5hmC, which could lead to erroneous conclusions regarding the relationship between different modified bases and their role in gene regulation. Therefore, we focused our discussion on studies that consider both 5mC and 5hmC. Several studies have investigated the DNA methylation landscape of certain subsets (Th1 and Th2) using Bisulfite sequencing and hydroxymethylation with immunoprecipitation techniques (35). The authors highlighted how Th1 and Th2 cells could be distinguished by 5mC and 5hmC changes in important identity genes. We furthered this work by mapping 5mC and 5hmC in important loci for these and additional Th cell subsets (Th17, Th17–no TGF-β and Tregs) and identified unique 5mCpG patterns for each subset and thus was able classify them by their epigenetic landscapes. Importantly, by characterizing methylation profiles in our identity gene panel, we were able to discriminate the unique Th17–no TGF-β cells, indicating that an epigenetic approach has the potential to identify additional pathogenic and intermediate Th cell subsets. This is in line with other studies showing that TGF-β in cultured cells may influence methylation profiles through cell division (36). Moreover, our own findings confirm that large differences in 5mCpG can also be identified in vivo in cells that escape TGF-β control. Emerging evidence suggests that metabolomics and epigenetics work in concert to alter Th cell identity (1).

DNA hydroxymethylation is an emerging marker of cellular identity and has been studied largely in the context of cellular differentiation (10, 37). For detection of multiple modified bases (5mC, 5hmC, etc.), most techniques require samples to be split, and different modified bases to be detected separately. In the current study, we took advantage of nanopore sequencing’s ability to determine the 5mCpG and 5hmCpG simultaneously. We found that 5hmCpG across a selection of genes was able to classify T cell subsets to a certain extent, and was largely concentrated in “active” genes. Our findings are in accordance with recent studies also highlighting that 5hmCpG in T cells is concentrated in the body of active genes (27, 35). Such studies indicated that terminally differentiated Th1 and Th2 cells have very little 5hmCpG (27) compared with CD4⁺, CD8⁺ or double-positive T cells (38). Our results obtained straight from in vivo differentiated gut Th17 cells suggest that the selected targets are able to discriminate between nonpathological and pathological Th17 cells. Of note, compared with the findings in in vitro–polarized T cells, the differences observed in vivo were stronger in number and magnitude. Not only were the differences in 5mCpG larger, but also we were able to observe higher basal levels of 5hmCpG and differences above 30% between conditions. For 5hmCpG in particular, our findings are in line with those suggesting that in vitro proliferation is associated with a passive loss of 5hmCpG during adaptation to culture conditions (39).

In addition to targeted DNA methylation analyses, we show here that off-target data can be used to effectively infer TF activity. Because T cell identity is defined by the synchronized activity of TFs, we illustrate how this can be exploited to complement locus-specific results without additional extraction or sequencing workload. We foresee that a combination of both on-target and off-target data will be a more robust marker of Th cell identity.

An important consideration for interpretation of the current study was that our samples were not purified after in vitro polarization and are known to be heterogeneous (cytokine treatment is not 100% efficient). However, we expect stronger discriminatory capacity in purified cell types. In addition, we selected important loci for Th cell subsets through literature review, which relied on known well annotated regions because there has yet to be a long read whole-genome sequencing study of T cell subsets. Although we used analysis of off-target loci to support our main findings and develop an unbiased profile, there are potentially additional genomic regions with better discriminatory capacity that are yet to be identified.

In our study, we used two different methods for targeted sequencing, CRISPR-Cas9 selection and UNCALLED adaptive sampling. The major difference between these two strategies is the moment where genomic selection occurs, either during library preparation in the case of CRISPR or during the sequencing process for UNCALLED. Although CRISPR-based libraries pass randomly through the nanopores in a given flow cell, adaptive sampling relies on real-time electronic selection of the molecules of interest. More precisely, UNCALLED converts stretches of raw signal into k-mers that are used as a query for an index search against the target database of genes of interest. These methods do not affect the calling of canonical and modified bases, and the downstream methods of analyses are indeed very similar. However, UNCALLED offers the advantage of a larger selection of targets (more than 100 in our hands) and a simplified library preparation. Moreover, other methods of adaptive sampling are currently available (40) or under development and will likely improve the coverage and efficiency of real-time sampling over presequencing target selection.

In conclusion, by mapping modified DNA bases (5mCpG and 5hmCpG) at single-molecule resolution, we show that specific epigenetic footprints are able to distinguish some of the major groups of in vitro polarized T cells (i.e., Th0, Th1, Th2, Tregs, and Th17s). Moreover, we show that intermediary phenotypes associated with dysfunctional immune responses can also be identified by targeted 5mCpG/5hmCpG profiling in vitro (Th17–no TGF-β cells) and in vivo (TGFbR-KO Th17 cells).

C.G. and H.H.-V. have received travel and accommodation support to attend conferences for Oxford Nanopore Technology. The other authors have no financial conflicts of interest.

We acknowledge the expert assistance of H. Tarayre V, S. Rodriguez, and E. Guillemot. We thank Profs. G. Stockinger and A.Y. Rudensky for providing IL-17–CRE and Foxp3-GFP mice, respectively. The help of the core facilities, the flow cytometry facility, and the mouse facility Plateforme Petit Animal du of the Cancer Research Center of Lyon is also acknowledged.