Early targeted DNA methylation profiling of CD04+/CD08+ T cells reveals pathogenic mechanisms in different stages of impaired glucose homeostasis Cardiovascular Diabetology

Background: Despite current intensive treatments, hyperglycemic patients have an unfavorable prognosis due to severe CHD and development of complications. The interplay between hyperglycemia and systemic inflammation can modify patterns of gene expression mainly affecting DNA methylation in promoter regions and, thus, endothelial damage. However, DNA methylome of CD04 + and CD08 + T cells, which play a relevant role in endothelial dysfunction has not been studied, especially during increasing hyperglycemia. Purpose: To identify differentially methylated regions (DRMs) in CD04 + and CD08 + T cells by comparing healthy subjects (HS) to Pre-Diab and type 2 (T2D) patients. This approach would investigate possible biomarkers useful to identify vascular damage already at Pre-Diab state. Methods: In this pilot study, we enrolled a subgroup of patients from our ongoing PIRAMIDE clinical trial (NCT03792607) including a total of 14 individuals classified in HS (n=2), Pre-Diab (n=6), and T2D (n=6). The DMRs were identified by using the reduced representation bisulfite sequencing platform (RRBS) which captures the majority of promoters and CpG islands. Big data analysis was performed by using the R package and ChromHMM algorithms. Results: Most of the total DMRs (30-35%) were in the promoter regions in increasing hyperglycemia vs HS. A global analysis of DMR-related genes overlapping between Pre-Diab and T2D patients showed a prevalence of DNA hypermethylation in both T cells. Interestingly, the secreted protein acidic and cysteine rich ( SPARC ) gene was annotated to the most hypomethylated-DMR in Pre-Diab and its methylation level gradually decreased in T2D patients. analysis and clinical data in different stages of impaired glucose homeostasis. Our pilot study established that the differential methylation of genes involved in T2D pathogenesis, such as SPARC , correlated with DBP, Creatinine, LVDd, LVSD, LAD, LVPWd, and AODd. Our data need to be confirmed in large multicenter in large multicenter studies. This evidence suggested a putative biomarker useful to early diagnose Pre-Diab patients and predict the presence/absence of CV and kidney complications.

may be a useful biomarker of vascular complications in Pre-Diab patients with a possible role for primary prevention. However, larger multicenter trials are needed to validate our epigenetic data in the clinical arena.

Background
According to the World Health Organization (WHO), the prevalence of prediabetes (Pre-Diab) metabolic state has increased in the past 10 years and its incidence will continue to increase in the coming years [1][2][3][4]. This intermediate state between normoglycemia and diabetes is associated with an increased risk of type 2 diabetes (T2D), cardiovascular diseases (CVDs) and dementia [1][2][3][4]. Currently, Pre-Diab affects about 40% of the adult population (1) and about 5-10% of people with prediabetes becomes diabetic every year, while the remaining part could present coronary heart disease (CHD) and CV complications [5,6]. Literature studies reported that 1) not all prediabetic patients progress towards T2D, 2) it is possible that no intervention is necessary in Pre-Diab state since there is no disease and 3) it could be associated with an increased risk of CVDs [1]. In order to obtain a better understanding of the mechanisms underlying progress of Pre-Diab toward T2D and CV complications, we studied some of the epigenetic modifications that can occur named DNA methylation [7]. Clinical evidence reported that peripheral insulin resistance, hyperglycemia and inflammation can significantly alter DNA methylation in early stages by providing a putative mechanistic link between glucose homeostasis imbalance, and vascular damage [8][9][10][11][12][13][14][15]. DNA methylation mainly occurs at cytosine bases located both in gene promoter and "CpG islands", which represent large regions with high frequency (about 50%) of cytosine-phosphate-guanine (CpG) dinucleotides [7,16,17]. Generally, methylation of CpG sites in gene promoter or associated CpG islands can modulate gene expression in a spatio-temporal manner [7,16,17]. We hypothesized DNA methylation changes can appear already at the Pre-Diab state or during the transition to T2D, providing putative useful early biomarkers for primary prevention.
Our group has a longstanding experience in epigenetics and CV diseases [18][19][20][21][22][23][24][25]. This pilot study is part of the clinical trial (NCT03792607) aimed at investigating early epigenetic changes in different stages of impaired glucose [26]. We performed a very complex DNA methylome analysis on both CD04 + and CD08 + T cells isolated from healthy subjects (HS), Pre-Diab and T2D patients (Figure 1). Our aim is to identify differentially methylated regions (DMRs) and annotated genes to clarify if these circulating cells may carry out detrimental signals underlying vascular damage. Indeed, as previously reported in acute coronary syndrome (ACS) patients [27], which result strongly associated with hyperglycemia (28), alterations of methylation signatures in both CD04 + and CD08 + T cells provided potential clinical biomarkers and therapeutic targets [29][30][31]. Since DNA methylation at regulatory regions is highly correlated with cell-specific patterns of repressive chromatin marks (7), we chose the reduced representation bisulfite sequencing (RRBS) platform which provides an enrichment of both promoters and CpG islands (32).
Thus, the identification of novel early and non-invasive molecular biomarkers by using liquid-based assays [33] could help physicians to select high risk hyperglycemic patients, allowing a selective phenotyping of the prediabetic patient and help to stratify the risk of developing CHD and vascular complications.

Patient enrolment
In this pilot study, we enrolled a subgroup of patients from our ongoing PIRAMIDE clinical trial (NCT03792607) (26) including a total of 14 subjects classified in HS (n=2), Pre-Diab (n=6), and T2D (n=6). Clinical characteristics were reported in Table 1. Our study population was recruited at Department of Advanced Medical and Surgical Sciences, University of Campania "Luigi Vanvitelli". Patients were clinically diagnosed as Pre-Diab and definite T2D based on clinical history, symptoms, and laboratory tests, in according to the current guidelines [34]. Pre-Diab was diagnosed by evidence of fasting plasma glucose of ≥5.6 mmol/L but <7.0 mmol/L (100-125 mg/dL; impaired fasting glucose, IFG), a 2-hour glucose of ≥7.8 mmol/L but <11.1 mmol/L during a 75 g oral glucose tolerance test (GTT) (140-199 mg/dL; impaired glucose tolerance, IGT), or a plasma hemoglobin (Hb) A1c of ≥5.7% but <6.4 %. T2D was diagnosed by evidence of IFG>7.0 mmol/L (>125 mg/dL), post prandial glycemia >11.1 mmol/L (>200 md/dL), and evidence of HbA1c>6.6%. Patients with known history of malignancy disorders, active infections, and chronic or immunemediated diseases were excluded from the study to avoid confounder effects. As controls, we selected subjects free of hyperglycemia and with no signs and symptoms of vascular damages or medication use. This study was approved by the local Ethical Committee (Protocol N. 114) and all patients signed a written informed consent. The study was conducted in accordance with the principles of the Declaration of Helsinki.

Clinical data analysis
Statistical analysis was performed to compare clinical characteristics for Pre-Diab and T2D groups, using R software (version 3.03). Continuous variables were expressed as mean ± standard deviation or standard error. Unpaired Student's t-test were used for comparison between two groups. Categorical variables were expressed as percentage and were compared using the Chi-Square test or the Fisher's exact test. A p-value of <0.05 was considered significant ( Table 1).

Cell isolation and DNA extraction
Peripheral venous blood samples (25mL) were collected from HS and patients and processed immediately after blood draw. Peripheral blood mononuclear cells (PBMCs) were separated from whole blood by Ficoll gradient using Histopaque®-1077 (Sigma-Aldrich), according to manufacturer's instructions. Starting from a total amount of 5x10 7 /mL of PBMCs, two leukocyte subsets, including CD04 + and CD08 + T cells, were separated and purified according to standard protocols for magnetic-activated cell sorting by using EasySep TM Human Naïve CD04 + T Cell Enrichment Kit (Stem Cell) and EasySep TM Human Naïve CD08 + T Cell Enrichment Kit (Stem Cell), respectively. Samples of genomic DNA were immediately isolated from isolated CD04 + and CD08 + T cellsfor each study participant by using DNeasy Blood & Tissue Kit (Qiagen) according to manufacturer's indications. DNA concentration and purity were determined by using NanoDrop spectrophotometer (Thermo Scientific) and evaluating both A260/A280 and A260/A230 standard absorbance ratio. Finally, DNA integrity was checked on 1% agarose gel (Supplementary Figure 1).

Library preparation
Sequencing was performed at the Genomix4Life S.r.l. (Salerno, Italy) [35]. Briefly, 2 μg of gDNA were used for each library preparation. DNA samples were digested by MspI restriction enzyme and purified with the GeneJet PCR Purification Kit (Thermo Fisher Scientific). All libraries were prepared by TruSeq Library Prep Kit (Illumina) and bisulfite conversion was obtained using the EZ DNA Methylation-Gold Kit (Zymo Research). DNA amplification reaction was performed using PfuTurboCxHotstart DNA Polymerase (Agilent Technologies, USA) and the amplified fragments were purified by AMPure XP Beads.
Finally, they were quantified by the Agilent 4200 TapeStation (Agilent Technologies). Each DNA library was analyzed by paired-end sequencing read (2 × 75 cycles) on Illumina Nextseq 500.

DNA sequence processing and alignment
Raw reads were assessed for quality by using FastQC (v011.8, Babraham Bioinformatics, UK) and trimmed to remove Illumina adaptors and low-quality reads using TrimGalore (v0.6.3, Babraham Bioinformatics, UK) with the default settings. A special setting designed for trimming RRBS data in Trim Galore was also used in order to trim the first 2 bases from the 3' end of the sequence, artificially added during RRBS library preparation. Trimmed reads were then mapped to the in silico bisulfite converted human reference genome (GRCh38) by using Bismark v.0.22.1 [36] with default parameters, the aligner Bowtie2 (v2.3.5.1) (37) and non directional parameter applied. Only the reads that were uniquely mapped were used to make cytosine methylation calls using bismark methylation extractor. On average across the 14 samples investigated for CD04 + T cells, 20037853 reads were obtained, 19878036 reads were retained after quality filtering, 10581587 reads were uniquely mapped to the reference genome with 53% mapping efficiency and the proportion of C methylated in a CpG context was about 38%. On average across the 14 samples investigated for CD08 + T cells, 20336957 reads were obtained, 20106684 reads were retained after quality filtering, 10964908 reads were uniquely mapped to the reference genome with 54% mapping efficiency and the proportion of C methylated in a CpG context was about 36%. Only methylation within CpG context was considered for further analysis.

Identification of differentially methylated regions (DMRs)
We identified DMRs by using the R package (v1.10.0) [38]. For both T cell types, we explored and compared three DNA methylation profiles: 1) HS vs Pre-Diab patients, 2) HS vs T2D patients and 3) Pre-Diab vs T2D patients. To prevent PCR bias and increase the power of the statistical test, CpG sites covering less than 10 reads or more than the 99.9 th percentile of coverage distribution in each sample were filtered out. Coverage values between samples were normalized as by default. Read coverage per base and correlation plots were calculated and displayed in Supplementary Figure 2. To define DMRs, we used a tiling window of 1 Kb. DMRs among the three groups were defined as regions with more than 25% methylation differences (|ΔM|) and q-value <0.01, after applying logistic regression by using the SLIM method to correct p-value for multiple hypothesis testing.
Positive and negative values indicate hyper-and hypomethylation in patients, respectively.

Gene annotation and functional analysis
By using the R package ChIPseeker (v1.20.0) [39], we represented the % of DMRs located into promoters, coding sequences, introns, distal intergenic regions, 5' UTR and 3' UTR in all three datasets, for both T cells (Supplementary Figure 3). In order to characterize the genes annotated to differentially methylated sites and regions, gene ontology (GO) and pathway analyses by using g:Profiler web server (http://biit.cs.ut.ee/gprofiler/index.cgi) were carried out [40]. This software calculates a minimum hypergeometric p-value for each nominated pathway, by placing greater weight on top-ranked genes where enrichment is the stronger. This strategy allows the detection of the most representative pathways of the entire gene list. For statistical significance, we considered all genes of the Human genome in the Ensembl database and applied the Bonferroni correction with a significant threshold <0.01 to adjust for multiple testing. GO analysis covered biological process (BP), molecular function (MF) domains and human phenotype (HP), whereas for biological pathways we considered both KEGG and Reactome databases.

Correlation analysis
After DMRs annotation, we considered methylation levels of relevant common DMR-related genes from both T cell types for correlation analysis. Significant association between their methylation status (hyper-and hypo-) and quantitative clinical parameters/laboratory tests for both T cell types and from each disease group was calculated by Spearman's Correlation, using as threshold corr> 0.5.

Chromatin state discovery and characterization
We used the ChromHMM (v1. 19) software (41) to characterize epigenomic regions according to 18 chromatin states, defined from both CD04 + and CD08 + T cells and grouped in HS, Pre-Diab and T2D. ChromHMM acquires chromatin state signatures by using a multivariate Hidden Markov Model (HMM) that clearly fits the combinatorial and spatial presence or absence of each chromatin state. From these signatures, ChromHMM generates a genome-wide annotation for each cell type and condition by calculating the most probable state in each genomic segment. Then, we annotated two methylated groups by applying the R package ChIPseeker and we found 564 genes annotated to hypo-states and 8389 annotated to hyper-states. Finally, we performed functional analysis for only promoters regions (<=1Kb), by using g:Profiler web server.
Raw data have been deposited in the NCBI Gene Expression Omnibus (GEO) database under the accession number PRJNA600866 (https://www.ncbi.nlm.nih.gov/sra).

Clinical characteristics of the study participants
Our study population was characterized by 57% (

Unique DMR trend in HS vs Pre-Diab and T2D
In order to characterize the most significant DMRs in distinct Pre-Diab and T2D groups, we firstly showed distribution values for each sample. Annotated DMRs were mapped according to their distance from established CpG islands. Notably, DNA methylation dynamics were rather concentrated among CpG islands located in the promoter regions  Table 2. Moreover, we identified the genomic locations most impacted by changes in DNA methylation between HS and increasing hyperglycemia.
Then, we discerned hypo-and hypermethylated DMRs-related genes (Supplementary Figure 5). In general, we observed that the number of hypermethylated DMRs was greater than those hypomethylated in CD04 + T cells; on the contrary, CD08 + T cells showed a higher number of hypomethylated DMRs than hypermethylated ones.

Common DMR trend in CD04 + T in HS vs Pre-Diab and T2D
We identified annotated genes common to Pre-Diab and T2D patients in a total number of 53, of which 11 were hypo-and 42 hypermethylated (Supplementary Figure 6).
Interestingly, from heatmap we noticed that there were some clusters of DMR-related genes, which retained the same hyper-or hypomethylation level in increasing hyperglycemia, as depicted in red boxes (Supplementary Figure 6). We focused on the top 18 highly significant DMR-related genes, which were shared in HS vs increasing hyperglycemia, of which 13 were hyper-and 5 were hypomethylated (Figure 2) (Table   3)

Common DMR trend in CD08 + T of HS vs Pre-Diab and T2D
We identified a total of 60 annotated genes, of which 40 were hypo-and 19 hypermethylated in HS vs hyperglycemia (Supplementary Figure 6). We focused on the top 20 highly significant DMR-related genes shared in HS vs increasing hyperglycemia, of which 7 were hyper-and 13 were hypomethylated (Table 4)

Analysis of CD04 + T cells
The top biological functions identified by g:Profiler and the number of genes involved are summarized in Supplementary Table 8. DMR-related genes of CD04 + T cells in Pre-Diab patients were mainly involved in the early damages to the micro-domains leading to abnormalities of eye morphology, physiology and movement as well as neurogenesis.
Moreover, a major involvement of nervous system differentiation was reported.

Analysis of CD08 + T cells
The top biological functions identified by g:Profiler and the number of genes involved are summarized in Supplementary Table 9. In Pre-Diab patients, DMRs induced abnormalities of central nervous system mainly involving high mental function and brain morphology as well as eye abnormalities. Moreover, it raised early damages in cardiac muscle tissue, development of blood vessels, as well as liver, limb and muscle physiology.

Analysis of CD04 + T cells
We evaluated the functional characteristics and signaling pathways associated to DRMrelated genes, which are shared by patients with increasing hyperglycemia (Supplementary Table 10). From results, the binding protein was the most prominent in the molecular process group, followed by catalytic activity (mainly hydrolase and transferase enzymes), nucleic acid binding, and transcription process. Interestingly, the human phenotype group highlighted abnormalities of vasculature, mainly aortic morphology and cardiac system. Moreover, abnormalities in eye, digestive system, muscle skeletal system (mainly hypotonia), immune system and liver were predicted.

Analysis of CD08 + T cells
A detailed GO analysis of significantly DRM-related genes in increasing hyperglycemia was summarized (Supplementary Table 11). At biological process level, a major regulation of immune system and neuron development was observed followed by hematopoiesis and response to lipid. KEEG highlighted the involvement of the insulin signaling pathway as well as REACTOME pointed to signal transduction (cytokines, interleukins, receptor tyrosine kinases), immune system, metabolism, transcription regulation and neural system, mainly PPIs at synapses.
We applied ChromHMM, a machine learning approach evaluating epigenomic information (called marks) across multiple cell types (CD04 + and CD08 + T cells) and multiple conditions (HS, Pre-Diab, T2D). As reported, this method allows to recognize chromatin states, by identifying their genomic occurrences [36]. The combination of multiple marks and the relative genomic annotation can be highly informative of distinct biological functions. Starting from 18 emission states (Supplementary Figure 7) ( Supplementary Table 12), we selected only genomic regions associated to a decreased (state 1) or increased (states 7, 9) methylation level vs HS. Then, we found 600 genomic regions associated to hypo-state 1 and 13947 associated to hyper-states 7 and 9 from Pre-Diab to T2D. Interestingly, we observed that the most states were distinctly associated with hypermethylation status during increasing hyperglycemia. As a consequence, we supported that a continuous glucose intake could induce a progressive DNA hypermethylation in both T cells.
Since hypermethylation of CpG islands clustered around transcription start sites (TSS) usually modulate gene expression [37], we considered only methylated regions annotated to promoters (<=1Kb). Then, we performed pathway analysis for these two groups of genes, by using web tool g:Profiler. We founded a total of 61 genes and 1198 genes annotated to hypo-state and hyper-states, respectively (Figure 3). We noted that genes associated to hypo-state enriched "phenotypic abnormalities". In particular, several genes, such as LRP5, ANGPTL6, SDHC and CCM2, characterized abnormalities of "CV system" and "blood circulation". Instead, for genes associated to hyper-states, we found an enrichment in many KEGG and REACTOME pathways and numerous "human phenotypes". In particular, KEGG analysis revealed three crucial signaling pathway, including MAPK, which plays an important role in the "regulation of CV diseases" and "cerebrovascular diseases" [42], PPAR, which is involved in "myocardial metabolism", [43] and ERBB, involved in "diabetes-induced cardiac dysfunction" [44]. From REACTOME, we found 54 genes, such as CD36, PPP2R5C and SOS1 involved in the pathogenesis of "insulin resistance", "glycemic control of T2D" and "CV diseases". Moreover, we noticed a lot of genes involved in "innate and adaptive immune system pathways". Finally, also through GO analysis, we found alterations in numerous "human phenotypes". A total of 51 genes characterized the increased "inflammatory response". In particular, some cytokine genes, such as interleukin 1 (IL-1), several complement cascade members, such as C5 and C4B showed a pathogenic role in "arterial hypertension", whereas SMAD3 mediated "diabetic cardiac hypertrophy". A large common set of genes was involved in "abnormality" of vascular (n=71), heart (n=45) and systemic arterial (n=28) morphology. Instead, a small group of these common genes, including SMAD3, NDE1, RAD51C, ZMPSTE24, FANCD2, PPARG, MYLK, GTF2I and FANCG, characterized also "abnormal carotid artery morphology".
Finally, another smaller group of genes, containing ABCC8, SAG and UCP2, was associated to "abnormal insulin level".

Discussion
Our study showed that: 1) the number of hypermethylated DMRs was greater than those hypomethylated in CD04 + T cells during increasing hyperglycemia, whereas the opposite trend was observed in CD08 + T cells; 2) functional analysis of unique DMR-related genes revealed that in CD04 + T cells early modifications of DNA methylation were already evident in Pre-Diab patients leading to a possible microvascular damage at eye level.
Moreover, we did not observe a predictive epigenetic role of CD04 + T cells in vascular damage for T2D patients. In contrast, we observed that different methylation profiles in CD08 + T cells were involved both in micro-and macrovascular damage from Pre-Diab to T2D patients (Figure 4); and 3) circulating T cells showed a set of overlapping DMRrelated genes with increasing or decreasing levels of DNA methylation from Pre-Diab to T2D patients vs HS. Prediabetes, as intermediate state between normo-and hyperglycemia, represents a high risk for T2D onset and vascular damage in asymptomatic patients [1][2][3][4]; however, there are no stringent diagnostic criteria to define when and what pharmacotherapy may aid to prevent CV dysfunction at individual level. Since inflammation strongly characterizes progression from HS to Pre-Diab, T2D and insulin therapy, additional biomarkers may improve primary prevention of CVDs.
The strengths of our study are that we evaluated and compared the differential DNA methylation profiles of two T lymphocyte subpopulations, CD04 + and CD08 + T, which may contribute to early vascular damage in different stages of impaired glucose homeostasis.
Recently, a DNA methylome analysis reported a predominant contribution from CD04+ and CD08+ T cells in regulating expression levels of IL6R, FASLG, and CCL18 genes in ACS patients vs HS suggesting a significant role in disease pathogenesis [26]. Since pathogenesis of ACS is related to vasculature damage and diabetes onset [27], we hypothesized that DNA methylation changes in both CD04 + and CD08 + T cells may reveal early molecular signals of CV dysfunction in Pre-Diab to T2D.
Previous DNA methylome analysis were focused on peripheral blood mononuclear cells (PBMNCs) [13] and tissue biopsies [14,15,45,46] isolated from T2D patients vs HS, without considering prediabetic state. Moreover, PBMNCs and tissue biopsies are characterized by cell heterogeneity which does not fit with the cell-specific DNA methylation patterns limiting the identification for starting sites of disease pathogenesis. Another strength of our RRBS analysis is that we specifically analyzed DMRs rather than methylation level of single CpG dinucleotides since DMRs can really control spatio-temporal gene expression.
Remarkably, selected DMRs have the most statistical power and by-pass putative effects of genetic polymorphisms during epigenome-wide association studies [47].
Our epigenetic trajectories suggested that hyperglycemia might early affect the interactome of both CD04 + and CD08 + T cells already in Pre-Diab patients by regulating DNA methylation at different set of genes. Interestingly, SPARC resulted the most hypomethylated gene in Pre-Diab and its methylation level gradually decreased in CD08 + of T2D patients. SPARC encodes for a multifunctional protein modulating the interaction between cells and the extra-cellular matrix (ECM) by the regulation of collagen and vitronectin [48,49]. SPARC is involved in many process including wound healing, inflammation, angiogenesis, remodeling of cardiac muscle and modulation of growth factor signaling [48,49]. Moreover, SPARC is expressed in adipocytes and pancreatic cells with profibrotic effects mainly linked to obesity and diabetic retinopathy and nephropathy [48].
Since the negative correlation between DNA methylation promoter and gene expression, we could expect gradually increased levels of SPARC proteins from Pre-Diab to T2D patients. Previous studies reported that increased plasma levels of SPARC protein were correlated to inflammation, insulin resistance and dyslipidemias in gestational diabetes [49]. Moreover, SPARC protein was increased in plasma of hyperglycemic patients, also correlating with early nephropathy in T2D [50]. This evidence may fit with our expected trend, for which DNA hypomethylation of SPARC promoter could increase levels of gene expression in increasing hyperglycemia. Targeted quantification of SPARC mRNA levels should be performed to confirm whether these changes of DNA methylation really lead to modulation of gene expression in CD08 T cells and/or other tissues isolated from Pre-Diab and T2D patients.

Conclusions
This is the first molecular-bioinformatic approach combining RRBS DNA methylome analysis and clinical data in different stages of impaired glucose homeostasis. Our pilot study established that the differential methylation of genes involved in T2D pathogenesis, such as SPARC, correlated with DBP, Creatinine, LVDd, LVSD, LAD, LVPWd, and AODd. Our data need to be confirmed in large multicenter in large multicenter studies. This evidence suggested a putative biomarker useful to early diagnose Pre-Diab patients and predict the presence/absence of CV and kidney complications.

Declarations
Ethics approval and consent to participate: University of Campanian "Luigi Vanvitelli" Ethical Committee (Protocol N. 114);

Cardiovascular Diabetology
Availability of data and materials: All data and materials are available following a specific request from the Editor;

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to