Comprehensive analysis of DNA methylation and gene expression in orally tolerized T cells

T cell anergy is known to be a crucial mechanism for various types of immune tolerance, including oral tolerance. The expression of several anergy-specific genes was reportedly up-regulated in anergic T cells, and played important roles in the cells. However, how the genes were up-regulated has not been understood. In this study, we comprehensively analyzed the altered gene expression and DNA methylation status in T cells tolerized by oral antigen in vivo. Our results showed that many genes were significantly up-regulated in the orally tolerized T cells, and most of the genes found in this study have not been reported previously as anergy related genes; for example, ribosomal protein L41 (FC = 3.54E06, p = 3.70E-09: Fisher's exact test; the same applies hereinafter) and CD52 (FC = 2.18E05, p = 3.44E-06). Furthermore, we showed that the DNA methylation statuses of many genes; for example, enoyl-coenzyme A delta isomerase 3 (FC = 3.62E-01, p = 3.01E-02) and leucine zipper protein 1 (FC = 4.80E-01, p = 3.25E-02), including the ones distinctly expressed in tolerized T cells; for example, latexin (FC = 3.85E03, p = 4.06E-02 for expression; FC = 7.75E-01, p = 4.13E-01 for DNA methylation) and small nuclear ribonucleoprotein polypeptide F (FC = 3.12E04, p = 4.46E-04 for expression; FC = 8.56E-01, p = 5.15E-01 for DNA methylation), changed during tolerization, suggesting that the distinct expression of some genes was epigenetically regulated in the tolerized T cells. This study would contribute to providing a novel clue to the fine understanding of the mechanism for T cell anergy and oral tolerance.


Introduction
Oral administration of food antigens is known to induce oral tolerance, and T cell anergy is reported as a major mechanism of oral tolerance as well as other various types of immunological tolerance [1][2][3]. Anergic T cells do not respond to the relevant antigen stimulation, while surviving for a long period of time. Although many studies have previously reported that the expression of several anergy-specific genes was up-regulated in anergic T cells [4][5][6][7], the mechanism for the regulation of their expression remains unknown.
As described above, the increased expression of anergy-specific genes is maintained over a long term [4][5][6][7]. Therefore, it has been suggested that some epigenetic regulations may be involved in the regulation of anergy-specific genes [8]; however, there is little evidence to a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 support this proposal. However, given that there are numerous genes showing altered expression levels in anergic T cells, it is unlikely that all the genes are independently and epigenetically regulated. Therefore, we are considering that only a few anergy-specific genes are epigenetically regulated and control the expression of other anergy-specific gene expressions. Indeed, in the case of other T cell subsets, a certain critical gene acts as a master regulator for each respective subset; for example, T-bet, GATA-3, RORγt and Foxp3 for Th1, Th2, Th17 and Treg cells [9][10][11], respectively. It is expected that the induction of T cell anergy is also regulated by a putative master regulator. In addition, some of the former four have been suggested to be epigenetically regulated [12], suggesting that epigenetic regulation is critical to controlling the regulators' expression.
We had performed a transcriptome analysis and a genome-wide DNA methylation analysis of T cells that were anergized in vitro using the next-generation sequencing technique [13]. Consequently, we found that the expressions of many genes were changed by anergy induction; for example, neuritin 1 (FC = 2.82, p = 1.08E-03: Fisher's exact test; the same applies hereinafter) and acid-sensing (proton-gated) ion channel 3 (FC = 2.72, p = 7.79E-07), and that the DNA methylation status of some of those genes was also changed; for example, neuritin 1 (FC = 3.00E-01). Based on the results of the study, we have not identified any master regulators of anergic T cells yet; however, the observations do indicate that the induction of T cell anergy is regulated by some epigenetic mechanisms.
In the present study, we performed a transcriptome analysis and a genome-wide DNA methylation analysis using T cells tolerized by oral tolerance as well as the previous study using in vitro anergized T cells. We considered that the orally tolerized T cell population included anergic T cells to a certain extent. In the current study, we carried out the study in two purposes; first, we aimed to identify candidates for the master regulator of anergic T cells induced by oral tolerance. Second, we aimed to confirm if the candidates would correspond to those obtained from the previous in vitro study [13]. Our present results provide several novel evidences about the features of in vivo tolerized T cells and some important clues to understand the mechanism for anergy induction.

Mice
Ovalbumin (OVA)-specific TCR-transgenic DO11.10 mice were obtained from the Jackson Laboratory (Bar Harbor, ME, USA). Their offspring were used at 5-7 weeks of age. The T cells of these mice recognize OVA 323-339 restricted to I-A d . The TCR-transgenic mice were bred in the animal facility of our university and were maintained on irradiated food and autoclaved distilled water. All mice were maintained and used in accordance with the guidelines for the care and use of experimental animals of Tokyo University of Agriculture and Technology. The procedures of this study were approved by the animal care and use committee in Tokyo University of Agriculture and Technology (26-8, April 1 st , 2014; 27-3, April 17 th , 2015; 28-3, April 19 th , 2016; 30-6, April 6 th , 2018). We checked the mice almost every day and confirmed no adverse clinical signs through the experimental periods. Splenocytes were prepared from the mice immediately after sacrificing them by cervical dislocation.

Induction of oral tolerance
Oral tolerance was induced in mice against OVA by feeding them 20% OVA-containing diet (Oriental Yeast Co., Ltd, Tokyo, Japan) for 7 days.

Restimulation
The splenocytes from the oral tolerance-induced mice were restimulated by incubation (2 × 10 6 cells/mL) in the presence of OVA (100 μg/mL) for 3 days. OVA was obtained from Sigma (St Louis, MO, USA).

ELISA
The supernatant was collected from each well 3 days after the restimulation procedure and was assayed by ELISA for its cytokine concentration. To measure IFN-γ and IL-4 concentrations, Maxisorp immunoplates were coated with purified R4-6A2 rat anti-mouse IFN-γ or purified 11B11 rat anti-mouse IL-4 antibody (BD Pharmingen, San Diego, CA, USA). Each sample and the standard were added after washing and blocking. The bound IFN-γ or IL-4 was detected using XMG1.2 biotinylated rat anti-mouse IFN-γ or BVD6-24G2 biotinylated rat anti-mouse IL-4 (BD Pharmingen) before incubating with alkaline phosphatase-streptavidin. The enzyme substrate-p-nitrophenol phosphate-was added, and the absorbance was determined at 405 nm.
To evaluate cell proliferation, splenocytes were also collected from each well 3 days after the restimulation procedure, and the cell proliferation level was measured by Cell proliferation ELISA, BrdU kit (Roche Diagnostics, Basel, Switzerland).

Purification of T cells
T cells were purified from cultured splenocytes (for transcriptome analysis) or fresh splenocytes (for DNA methylation analysis and RT-PCR) by the MACS separation system (Miltenyi Biotec, Bergisch Gladbach, Germany) using CD90.2 micro beads and an LS column.

Transcriptome analysis
To prepare cDNA library, total RNA was isolated from normal or tolerized T cells 10 h after restimulation using the Tripure isolation reagent (Roche). To avoid DNA contamination, total RNA was treated with DNase (Takara Bio, Kusatsu, Japan). The cDNA libraries were synthesized using mRNA template and random hexamer primers; then, a custom second-strand synthesis buffer (Illumina, San Diego, CA, USA), dNTPs, RNase H, and DNA polymerase I were added to initiate the second-strand synthesis. After a series of terminal repair, a ligation and sequencing adaptor ligation, the double-stranded cDNA libraries were completed through size selection and PCR enrichment. The qualified libraries were fed into an Illumina HiSeq sequencer after pooling according to its effective concentration and the expected data volume.
The quality of the sequencing reads was checked by FastQC (Babraham Bioinformatics, Cambridgeshire, UK), and any low-quality sequences were trimmed by a custom-made script with Perl, v5.10.1. The trimmed sequence reads were subjected to further analyses. The reads for RNA-Seq were mapped onto mouse cDNA reference sequences (Mus_musculus.GRCm38. cdna.abinitio.fa; ftp.ensembl.org) by Bowtie2 (version 2.0.6) [14]. The mapped reads were counted and calculated as FPKM. All statistical analyses were performed by GNU bash (version 4.3.48(1)-release) and R (version 3.4.3). Statistical significance of differential expression was represented by p-value with Fisher's exact test. Biological annotation and analyses, including Gene Ontology (GO) analysis and KEGG pathways analysis, were performed with DAVID Bioinformatics Resources 6.8 software (https://david.ncifcrf.gov/).

RT-PCR
cDNA was synthesized from total RNA obtained from restimulated normal or tolerized T cells by using a Transcriptor First Strand cDNA Synthesis kit (Roche), and qPCR was then performed (Thermal Cycler Dice Real Time System II, Takara Bio, Kusatsu, Japan) by using KAPA SYBR FAST qPCR kit (KAPA Biosystems, Wilmington, MA, USA). The level of expression of GRAIL was standardized by that of β-actin. The specific primers for each gene are listed as follows: β-actin, 5'-GTGGGCCGCTCTAGGCACCAA-3' and 5'-CTCTTTGATGTCACG CACGATTTC-3'; GRAIL, 5'-AGAGAGAGGGGCTTCTGGAG-3' and 5'-CGATGACCATT GTGACTTGG-3'. The statistical analysis was performed by Student's t test.

Bisulfite treatment of genomic DNA
To analyze the changes in the DNA methylation status owing to the induction of oral tolerance, genomic DNA was prepared from the normal and tolerized T cells by the phenolchloroform method. The DNA sample was treated with bisulfite using a MethylEasy Xceed Rapid DNA Bisulphite Modification Kit (Takara Bio). The completion of the bisulfite treatment was confirmed by PCR using Epitaq HS (for bisulphite-treated DNA, Takara Bio) according to the manufacturer's protocol.

DNA methylation analysis
The DNA samples that were or were not treated with bisulfite were sequenced by an Illumina GAIIx sequencer using an TruSeq DNA Methylation kit (Illumina) according to the manufacturer's protocol.
Sequence reading of bisulfite-genome sequencing for DNA methylation analysis was performed with mouse genomic DNA reference sequences (Mus_musculus.GRCm38.dna.toplevel. fa; ftp.ensembl.org) by Bismark (Babraham Bioinformatics, Cambridgeshire, UK) [15]. Annotation and statistical analyses for DNA methylation analysis were performed by using a gff file (Mus_musculus.GRCm38.84.gff3; ftp.ensembl.org) by custom-made scripts with GNU bash (version 4.3.48(1)-release), Perl (v5.22.1), and R (version 3.4.3). Statistical significance of differential methylation level was represented by p-value with Fisher's exact test. Biological annotation and analyses, including Gene Ontology (GO) analysis and KEGG pathway analysis, were performed by DAVID Bioinformatics Resources 6.8 software. The methylation rate of cytosines within 10 kb upstream from the transcription start site (TSS) of each gene was analyzed.

Orally tolerized T cells specifically expressed various genes
To identify the genes specifically expressed in orally tolerized T cells, the cells were subjected to transcriptome analysis. We first confirmed the induction of tolerization in the splenocytes of the orally tolerized mice [1,4]. Our results showed that all four parameters-proliferation (Fig 1A), IFN-γ ( Fig 1B) and IL-4 ( Fig 1C) production were significantly suppressed, and GRAIL expression (Fig 1D) was significantly up-regulated-were changed in the tolerized T cells. Next, we prepared the mRNA samples and libraries for transcriptome analysis using T cells purified from the antigen-restimulated splenocytes. We confirmed that the purity of the T cells was approximately 90%. We also confirmed that the frequency of OVA-specific T cells in DO11.10 was approximately 65%. We analyzed three samples obtained from different mice for each group and showed the data regarding the number of reads obtained from DNA sequencing for transcriptome analysis (S1 Table). We defined the altered expression of genes caused by tolerization according to the following criteria: increased by 2.0-fold or more, or decreased to 0.5-fold or less. Our results showed that we could detect the expression of 32926 genes, and that the expressions of 13.02% of the genes were changed in orally tolerized T cells. Among them, 2034 genes; S2 Table (6.18% of them) were increased and 2286 genes (6.94%) were decreased (Fig 2). To understand the biological features of orally tolerized T cells, we analyzed the functions of the 2034 genes that showed increased expression in the tolerized T cells using DAVID Bioinformatics Resources 6.8. Consequently, we found 858 genes on the DAVID database and categorized them into 107 clusters. We have shown the top 5 of them in Table 1. Although DAVID does not give us the name of each cluster, the results would show that clusters 1 and 2 included genes involved in regulation of gene expression by remodeling the chromatin structure, cluster 3 included genes involved in host defense against infections, cluster 4 included genes involved in general immune responses, and cluster 5 included genes involved in binding the metal ions. To identify enriched pathways, we performed the KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis using DAVID Bioinformatics Resources 6.8. We confirmed 1360 genes were assigned by the KEGG analysis among genes with up-regulated expression in the tolerized T cells by 2.0-fold or more. The result showed that several pathways including those involved in immune responses, for example 'Cytokine-cytokine receptor interaction' (Fold enrichment = 2.16, p = 1.14E-4), 'Chemokine signaling pathway' (Fold enrichment = 2.05, p = 1.79E-3) and 'TGF-beta signaling pathway' (Fold enrichment = 2.46, p = 8.75E-3), were changed in the tolerized T cells (Table 2). Next, we analyzed the genes upregulated in the tolerized T cells in higher detail. In Table 3, we show the genes selected among 121 genes with significance by Fisher's exact test by following 2 criteria: 1) the top 10 for fold change of expression, or 2) detected in all three tolerized samples. In addition to those genes, we also focused on transcription factors. This list includes several interesting genes, for exam-  Percentage of genes with up-or down-regulated expression during tolerization. Splenocytes from DO11.10 mice that had been fed OVAcontaining diet for 7 days were restimulated with OVA. Total RNA was prepared from the cells 10 h after the restimulation, and underwent transcriptome analysis. Genes with fold change rate �2.0 were considered to be up-regulated, while those with fold change rate �0.5 were considered to be down-regulated. The functions of genes with increased expression in the tolerized T cells were analyzed by DAVID Bioinformatics Resources 6.8. We found 858 genes on the DAVID database and categorized them into 107 clusters. We have shown the top 5 among them in this

The methylation status of various genes altered by tolerization
Then, we comprehensively analyzed the genome-wide methylation status in both control and orally tolerized T cells. We subjected the DNA samples prepared from the tolerized T cells to DNA sequencing after treating the DNA with bisulfite. We analyzed three samples obtained from different mice for each group and showed the data regarding the number of reads obtained from the sequencing in S3 Table. Tolerization-induced alteration in the methylation status of the genes was defined according to following criteria: increased by 2.0-fold or more, or decreased to 0.5-fold or less. Our results showed that among the identified whole genes (30136 genes), 4.13% (1245 genes) of the genes showed altered DNA methylation status, in orally tolerized T cells; among them, 3.18% of the genes (960 genes) showed increased and 0.95% (285 genes; S4 Table in red characters) showed decreased DNA methylation (Fig 3). Then, we used DAVID to analyze the functions of 285 genes with decreased methylation status in the tolerized T cells. We found 167 genes on the DAVID database and categorized them into 24 clusters. We have shown top 5 of them in Table 4. Clusters 1 and 2 included genes involved in regulation of gene transcription, cluster 3 included genes involved in reduction-oxidation reactions in the mitochondrion, cluster 4 included genes involved in methyltranferase activity, and cluster 5 included genes involved in cilium assembly. To identify enriched pathways, we performed the KEGG analysis. We confirmed 2317 genes were assigned by the KEGG analysis among genes with decreased DNA methylation in the tolerized T cells by 0.9-fold or less, using DAVID Bioinformatics Resources 6.8. The result showed that several pathways were changed in the tolerized T cells (Table 5). There are several pathways associated with diseases relevant to immune dysfunction, for example 'Maturity onset diabetes of the young' (Fold enrichment = 2.76, p = 2.07E-2) and 'Autoimmune thyroid disease' (Fold enrichment = 1.84, p = 3.69E-2). Next, we analyzed the genes with decreased methylation status in the tolerized T cells in more detail. We have shown the genes in Table 6. They were selected based on the following 4 criteria: 1) the fold change of the methylation rate of the gene decreased by 0.5-fold or less, 2) evaluated as significant by Fisher's exact test

Several genes up-regulated in orally tolerized T cells had a decreased DNA methylation status
Finally, we showed 12 genes in Table 7 as candidates for the master regulator of orally tolerized T cells, selecting according to following 4 criteria: 1) the fold change of expression increased by 2.0 or more (S2 Table), 2) the fold change of expression was evaluated as significant by Fisher's exact test (S2 Table), 3) the fold change of the methylation rate of the gene decreased by 0.9-fold or less (S4 Table), and 4) the fold change of the count of methylated cytosine either decreased or remained unchanged (S4 Table). The genes include Med9 (FC = 6.17E03, p = 1.59E-02 for expression; FC = 7.14E-01, p = 2.42E-01 for DNA methylation) involved in the regulation of transcription [16], Hist1h2ai (FC = 7.88E04, p = 4.75E-05 for expression; FC = 7.05E-01, p = 3.25E-01 for DNA methylation) involved in chromatin remodeling [17], and Lxn (FC = 3.85E03, p = 4.06E-02 for expression; FC = 7.05E-01, p = 4.13E-01 for DNA methylation) associated with inflammation [18]. On the other hand, most of them have not been reported yet if they have specific functions in anergic T cells.

Most cytosine in CpG islands were methylated
We calculated the rate of methylated cytosine in CpG, CHG, and CHH sequences. Our data showed that most of methylated cytosines were found in CpG islands and the methylation rate of CpG islands was approximately 75% in both normal and tolerized T cells (Table 8).

Discussion
T cell anergy is known to be an important phenomenon that is involved in several types of immune tolerance, such as self-tolerance, superantigen-induced tolerance, and oral tolerance [1][2][3]19]. Therefore, many studies have been conducted to clarify the molecular mechanisms for unresponsiveness of anergic T cells [20,21] albeit with little success. In the present study,  we aimed to identify the genes which are specifically expressed via epigenetic regulation in anergic T cells obtained from orally tolerized mice. Given that the expression of so many genes in T cells gets altered during tolerization, it is believed that one or few key molecule(s) get specifically expressed in anergic T cells, similarly with T-bet, GATA-3, RORγt, and Foxp3 for Th1, Th2, Th17, and Treg cells [9][10][11], respectively, and then regulate the expression of other genes. Some of these master regulators have been demonstrated to be controlled under epigenetic modification of DNA [12]. A comprehensive analysis of anergy-specific gene regulations would contribute to the fine understanding of molecular mechanisms of T cell tolerization. By transcriptome analysis, we detected 2034 up-regulated genes and 2286 down-regulated genes in orally tolerized T cells. This result shows that up-regulated genes were fewer than down-regulated genes, supporting the fact that various types of responses are suppressed in tolerized T cells. A similar result was observed in our previous study in which we had used T cells anergized in vitro [13].
The result of DAVID analysis indicated that clusters 1 and 2 included genes involved in chromatin remodeling. For example, PAPA-1 (listed in Table 3) is a subunit of INO80, which is known to induce chromatin remodeling by DNA binding [22]. This result suggests that the epigenetic regulation of the gene expression may be important to maintain the hyporesponsiveness of the cells. Interestingly, clusters 3 and 4 included various genes involved in positive  immune responses, such as cytokine and chemokine activities, suggesting that some T cells included in the tolerized T cell population are not completely unresponsive and can influence other cells around them in some specified way. This result is supported by the study of Zheng et al. showing that some anergic T cells could enhance peripheral tolerance via interactions with other immune cells [23]. In addition, some pathways were indicated by KEGG analysis to be altered in tolerized T cells, including ones important for immune regulation (Table 2). Therefore, dysregulation of such critical pathways would result in deficiency of T cell function. Particularly, TGF-β is known to be one of the important cytokine involved in oral tolerance. The change in TGF-β-related signaling pathway would play a key role in the tolerized T cells. All things considered, the results of DAVID and KEGG analysis suggests that the data obtained from the transcriptome analysis are reliable because the genes found to be up-regulated are valid as anergy-related genes. The genes listed in Table 3 included several interesting genes, such as CD52, Hspe1, Ire3, and C1qbp. These genes have been reported as negative regulators of immune responses. For example, the proliferation and IFN-γ production of T cells that highly express CD52 were shown to be inhibited [24]. The differentiation of T cells into Th17 cells and the production of IL-17 were reportedly enhanced in Ire3-deficient mice [25]. These evidences suggest that those genes contribute to maintaining the hyporesponsiveness of tolerized T cells. Table 3 also shows that several genes involved in transcription regulation were up-regulated in the tolerized T cells. Some of them (Rest and Lbh) directly regulate gene transcription as transcription factors, while others indirectly control it via binding to transcription factors (Yaf2, Chchd2, C1qbp, and Hexim1) or regulating the activity of RNA polymerase II (Med10, Med9, and Cks2). These genes may include the key molecules that regulate the expression of various genes specific to anergy induction in T cells.
The next-generation sequencing technology made it possible to comprehensively understand the DNA methylation status of animal cells. We detected about 20,000,000 methylcytosine residues, and most of them were found in CpG islands ( Table 8). The methylation rate of CpG islands was approximately 75%. This result shows that the data are reliable because it has been generally known to be 60%-90% in mammalian cells.
Our result showed that the genes with decreased methylation rate caused by oral tolerance were fewer than those with increased methylation rate (Fig 3), thus supporting the overall result demonstrating that genes with increased expression were fewer than those with decreased expression (Fig 2). As described above, these changes may be partly involved in hyporesponsiveness of tolerized T cells. For example, in a previous study, expression of IL-2 was reported to be regulated by the methylation status of specific CpG sites located in an area within approximately 600 bp of the promoter region [26].
Conversely, our result also showed that the methylation rate of 285 genes was down-regulated during tolerization, suggesting that the expression of some of these genes was up-regulated by epigenetic regulation. It is suggested that several specific methylation sites would particularly contribute to control the expression of genes. For example, the differentiation into Treg cell was critically determined by the demethylation of TSDR on the Foxp3 gene [12]. It is considered that the methylation status of CpG sites on the binding regions of RNA polymerases or transcription factors was particularly important for regulating gene expression. Therefore, in future studies, we need to analyze the change in methylation rate focusing on each cytosine to better understand the regulatory mechanism of the altered gene expression, using by other more efficient strategies, such as microarray-based DNA methylation analysis [27] and Agilent SureSelect Mouse Methyl-Seq [28]. We do not have sufficient capacity to perform that using the data of the genome wide DNA methylation analysis with bisulfite sequencing at present because such analysis needs a computer with higher performance and deeper coverage of data. On the other hand, the results obtained from this study according to the average change of methylation status within 10 kb upstream from the TSS are sufficiently reliable to achieve the aim of this study. It has been demonstrated that CpG islands at these regions have critical roles for regulating gene expression.
The result of DAVID analysis for demethylated genes indicated that clusters 1 and 2 included the genes involved in transcription regulation, such as Rfx3, Sp8, and Tfcp2l1, which are listed in Table 6. There is no previous report on the immunoregulatory functions of these genes, although Rfx reportedly works in pancreatic β cells [29] and Sp8 plays a role in the development of neural cells under epigenetic regulation [30]. On the other hand, Table 6 also includes several genes involved in immune responses. Sash3 regulates cytokine production [31], Gimap1 controls B and T cell differentiation [32], and Ptgis is associated with the suppression of NFκB activity [33]. Conversely, the functions remain unknown for approximately a half of genes listed in Table 6. The identification of their functions would greatly contribute understanding the molecular mechanism of T cell anergy.
The KEGG analysis indicated that DNA methylation status of genes involved in several pathways was altered in the tolerized T cells ( Table 5). Several of them suggested dysfunction of immune regulation, for example 'Maturity onset diabetes of the young' and 'Autoimmune thyroid disease'. These diseases are induced by abrogation of immune tolerance. Therefore, the alteration of the pathways in the orally tolerized T cells suggests the importance of those pathways in oral tolerance as well. However, how the change contributes to unresponsiveness of the tolerized T cells is needed to clarify in future investigation. Interestingly, the results of KEGG analysis shown in Tables 2 and 5 included a common pathway 'Transcriptional misregulation in cancer'. Supplementary figure (S1 Fig) illustrated the pathway. The genes mapped in the pathway were indicated as red stars. This pathway gives the cancer cells specific characteristics such as differentiation resistance and survival. These features are also likely to be seen in tolerized T cells.
We adopted several other criteria to pick up genes which were the candidates for specifically expressed genes via epigenetic regulation in orally tolerized T cells (Table 7) because we could not detect the expression of most genes listed in Table 6 by the transcriptome analysis. Briefly, we expanded target genes with the fold change of methylation rate decreased by 0.9-fold or less because only few specific methylation sites might be important to control the expression of the gene. We summarized the procedure to pick up the genes shown in Table 5 (Fig 4). These genes are expected to encode key molecules for maintaining the anergic state. Med9 was the only gene known as a transcription factor among the genes listed in Table 7. Med9 is a subunit of a mediator complex which interacts with the carboxyl domain of the major subunit of RNA polymerase II and regulates the expression of specific genes by acting as a transcription co-activator [16]. This evidence suggests the involvement of Med9 in the regulation of genes important for anergy induction. The experiments for revealing the function of Med9 in anergic T cells are ongoing. In addition, an experiment for obtaining more data regarding the expression of the genes listed in Table 6 is also now in progress. It is expected that we would find other candidate genes for key molecules through the experiment. Then, the experiment to analyze the function of the newly identified candidates as well as the genes in Table 7 will be performed.
One of the aims of this study was to compare the result from T cells tolerized in vivo with that from T cells anergized in vitro. As shown in Table 9, we found only few genes of which expression were up-regulated in both in vitro anergized and in vivo tolerized T cells. The number of the genes is less than expected, suggesting that the mechanism for inducing anergy would not always be the same for anergic T cells induced in vitro and in vivo. On the other hand, it cannot be denied that these few genes play critical roles in the induction of T cell anergy. In addition, the methylation status of Snrpf, Ccdc124, Med9, Gng5, and Gm15429  Data for in vitro anergized T cells were quoted from our previous study [13]. https://doi.org/10.1371/journal.pone.0229042.t009 among the genes listed in Table 7 was also comparably decreased in in vitro anergized T cells (Snrpf; FC = 0.655, Ccdc124; FC = 0.627, Med9; FC = 0.733, Gng5; FC = 0.333, Gm15429; FC = 0.464), suggesting that the change of methylation status of those genes would be regulated by a common mechanism in both cells, and such mechanism is crucial for inducing T cell anergy. Further research is needed to identify the common mechanism and the molecules involved in it. Several anergy-specific genes, such as GRAIL, Itch, Cbl-b, Egr2, and Egr3, reportedly play important roles in anergy induction [4-7, 34, 35]. Indeed, we confirmed the expression of GRAIL was certainly up-regulated by inducing oral tolerance (Fig 1D). In this study, we could analyze the methylation status of these genes, except for that of GRAIL; however, the methylation status of these genes did not change even after tolerization (Itch; FC = 1.441, Cbl-b; FC = 1.169, Egr2; FC = 1.261, Egr3; FC = 1.168), indicating that the expression of these genes is probably regulated by mechanisms other than DNA demethylation. Similar results were obtained in our previous study using T cells anergized in vitro as well. It is suggested that these genes are not the master regulators of anergic T cells, and other key molecules with epigenetically enhanced expression may control those genes.
We were wondering if the transcriptome and the comprehensive DNA methylation analysis could be applied for T cells tolerized in vivo because each T cell was tolerized in various degree in vivo compared to T cell anergized in vitro. In addition, T cells orally tolerized in vivo would include not only anergic T cells but also other distinct populations such as Treg cells. Contamination of those cells might have interfered the analysis even if they are not major population. Furthermore, purity of antigen-specific T cells used in this study was also a big concern. Contamination of antigen-nonspecific T cells and non-T cells might also have interfered the analysis. However, in this study, we could identify several genes specifically expressed and epigenetically regulated in orally tolerized T cells, suggesting that the methods used in this study was enough to find candidates for the master regulator while can be improved to obtain more precise data. Most of those genes identified in this study have not been previously reported as candidates for the master regulator of anergic T cells. These results could be obtained thanks to the comprehensive analysis using next-generation DNA sequencing technology. This strategy can be applied for the investigation of master regulators of other unique cells. This study would contribute not only to clarifying the mechanism of T cell anergy but also for understanding various biological events.
Supporting information S1 Table. The number of reads obtained from DNA sequencing for transcriptome analysis. Each cDNA library was prepared from the DNA samples obtained from the control and tolerized T cells. They were analyzed with Illumina HiSeq sequencer.