DNA methylation analysis of phenotype specific stratified Indian population

Abstract

Background

DNA methylation and its perturbations are an established attribute to a wide spectrum of phenotypic variations and disease conditions. Indian traditional system practices personalized medicine through indigenous concept of distinctly descriptive physiological, psychological and anatomical features known as prakriti. Here we attempted to establish DNA methylation differences in these three prakriti phenotypes.

Methods

Following structured and objective measurement of 3416 subjects, whole blood DNA of 147 healthy male individuals belonging to defined prakriti (Vata, Pitta and Kapha) between the age group of 20-30years were subjected to methylated DNA immunoprecipitation (MeDIP) and microarray analysis. After data analysis, prakriti specific signatures were validated through bisulfite DNA sequencing.

Results

Differentially methylated regions in CpG islands and shores were significantly enriched in promoters/UTRs and gene body regions. Phenotypes characterized by higher metabolism (Pitta prakriti) in individuals showed distinct promoter (34) and gene body methylation (204), followed by Vata prakriti which correlates to motion showed DNA methylation in 52 promoters and 139 CpG islands and finally individuals with structural attributes (Kapha prakriti) with 23 and 19 promoters and CpG islands respectively. Bisulfite DNA sequencing of prakriti specific multiple CpG sites in promoters and 5′-UTR such as; LHX1 (Vata prakriti), SOX11 (Pitta prakriti) and CDH22 (Kapha prakriti) were validated. Kapha prakriti specific CDH22 5′-UTR CpG methylation was also found to be associated with higher body mass index (BMI).

Conclusion

Differential DNA methylation signatures in three distinct prakriti phenotypes demonstrate the epigenetic basis of Indian traditional human classification which may have relevance to personalized medicine.

Introduction

DNA methylation at cytosine bases of CpG dinucleotide is an epigenetic phenomenon widely used to regulate gene expression in humans [1]. CpG islands are often present in 5′-upstream region of genes and methylation of these is associated with condensed chromatin and repressed gene transcription [2,3]. As a heritable change, DNA methylation plays an important role in embryonic development, genomic imprinting and X-chromosome inactivation [4,5]. Alterations in DNA methylation is known to be affected by number of environmental factors to influence human diseases and is also associated with genetic events that leads to several types of cancer [6]. Growing evidences suggest that enzymatic machineries such as DNA methyltransferases (DNMTs) are responsible for creating and maintaining patterns of DNA methylation in response to environmental factors and during developmental programs [7,8].

The broad classification of phenotypic variations based on physiological, anatomical and psychological evaluation is being practiced in traditional Indian system of medicine Ayurveda [9,10]. Ayurveda suggests that principles of these evaluations are the basis of individuals life course and health events established at birth or prakriti which can be identified as Vata, Pitta and Kapha phenotypes with discrete characteristics [11,12]. The human classifications based on the body constitution as Vata, Pitta and Kapha prakriti in Ayurveda has formed the basis of disease management and for practicing traditional personalized medicine [13]. Several studies have attempted to associate single nucleotide polymorphisms (SNPs) to identify genetic basis of prakriti classification such as for HLA alleles [14], fast and slow metabolizing activities of defined CYP2C19 SNPs [15], inflammatory and oxidative stress SNP markers in rheumatoid arthritis[16], SNPs for external environmental adaptation such as EGLN1 [17], risk factors for coronary artery diseases [18] and variations in platelet aggregation among prakriti phenotypes [19]. Gene expression analyses performed among the three prakriti have also identified biochemical and hematological differences [20]. Previously we have showed differences in expression of cluster differentiation markers CD14, CD25 and CD56 among prakriti phenotypes by immunophenotyping [21].

Phenotypic variations observed across vast human diversity cannot be explained by genetic factors alone and may encompass several epigenetic phenomenon. In a genome wide study, Heyn et al. [22] have successfully reported the identification of DNA methylation markers in three human populations namely, Caucasian-Americans, African-Americans and Hans-Chinese Americans and have suggested to contribute to natural human variations. Therefore, DNA methylation may be involved in human developmental plasticity, adaptation and manifestation of distinct phenotypes. We studied the association of natural DNA methylation variations in the context of human classification based on principles of prakriti to understand phenotypically described sub-groups. We report that prakriti based classification of individuals provide genome-wide distinct DNA methylation signatures which may help in deciphering the basis for human variation description.

Methods

Study participants

The protocol was approved by institutional ethics committee and samples were screened after obtaining written informed consent, from all the participants. The Institute of Ayurveda and Integrative Medicine (I-AIM), Bangalore, received a written ethical approval from the IAIM- Institutional Ethical Committee (IEC), Shri Dharmasthala Manjunatheshwara College of Ayurveda (SDMCA), Udupi, obtained written ethical consent from Shri Dharmasthala Manjunatheshwara College of Ayurveda, Udupi, Institutional Ethical committee and Sinhgad College of Engineering, Pune, obtained the ethical consent from the Ethics Committee of Interdisciplinary School of Health Sciences, Pune University.

Screening and selection for dominant prakriti phenotyping was performed by double blinded method in three different cross-sectional cohorts using classical ayurvedic parameters [9,10] and a validated software called Ayusoft (Ayusoft, C-DAC. Pune). Qualitative subjective assessment of prakriti by ayurvedic physician and was arbitrated by the use of Ayusoft, thus providing a quantitative approach of prakriti determination. Detailed background, phenotypes of prakriti and methodology for the selection of subjects are described in our previous publication (Additional file 1) [23]. Briefly, 3416 healthy male volunteers of age group between 20–30 years were screened after receiving the institutional consent. Among them, 971 volunteers showed one of the three prakriti predominantly (≥ 60%), from whom whole blood genomic DNA was isolated by standard procedure.

Global methylation analysis by reverse-phase HPLC

Reverse-phase high performance liquid chromatography (RP-HPLC) method was employed with minor changes for global methyl cytosine (mC) estimation as described previously [24,25]. In brief, 1.0 μg of genomic DNA was treated with 1U of DNaseI (New England Biolabs, USA), denatured at 100^°C for 10 minutes and rapidly cooled in ice followed by the addition of 1.0 U of Nuclease P1 (Sigma-Aldrich, Canada) and 2.0 U of calf intestinal phosphate (New England Biolabs, USA). The 5-mC content was estimated in duplicate by injecting the sample to the RP-HPLC using C18 (Grace Vydac, Hesperia, CA, USA) columns. Isocratic delivery of mobile phase consisted of 50 mM potassium dihydrogen phosphate (pH 3.5) and methanol mixed in 9:1 ratio respectively, with flow rate of 1 mL/min. The percentage of 5-mC was calculated by using the formula [(5-mC peak area)/(C peak area /5-mC peak area)] /100 (where C is Cytosine and 5-mC is 5-methyl Cytosine)].

Methylated DNA immunoprecipitation (MeDIP) Microarray

Briefly, 10 μg of whole blood genomic DNA was sonicated by using a 2 mm probe at amplitude 40 for 30 cycles with 15 seconds on and off, to obtain fragment sizes of 100–800 base pairs. About 3 μg of sheared DNA (INPUT) was immunoprecipitated overnight with anti-5-mathyl cytosine antibody (Diagenode, Belgium) along with positive and negative control DNA provided in the Diagenode kit. Enriched methylated fragments were purified by standard phenol-chloroform-ethanol precipitation [26]. The quality analysis of enriched immunoprecipitated (IP) methylated fragments was evaluated by real time PCR with positive, negative, GAPDH and AX1 primers (Diagenode, Belgium) [27]. INPUT and IP were differentially labeled with cyanine3 and cyanine5 (Amersham Biosciences, USA) by indirect method using Bio-prime Array CGH kit (Invitrogen, USA) [28,29]. Equal concentrations of labeled INPUT and IP DNA were co-hybridized onto the Agilent 244K Human CpG Island, high density oligonucleotide array at 65°C for 40 hours with the continuous rotation at 18 rpm as per the Agilent MeDIP protocol version 1.1 (Agilent Technology, USA). The slide was washed using wash buffers (Agilent Technologies, Santa Clara, CA, USA), dried and scanned using G2505B DNA microarray scanner (Agilent Technology, USA) with Sure Scan High resolution technology (Additional file 1).

Microarray Data Analysis

Total of 147 randomly selected male volunteers of dominant prakriti (Kapha = 52, Pitta = 48 and Vata = 47) were used for data analysis. Feature extraction software v10.1 was used to obtain background corrected, loess normalized, and logarithmic converted intensity values for green and red channels for each probe. Logarithmic difference of red to the green signal is the measure of methylation, generally represented as log ratio or (INPUT/IP) or “M-value” (methylation) for a given probe. Subsequent quality analysis was performed using R v2.15 with Bioconductor package limma. Quantile normalization was adopted for between array correction to remove batch effect and to equalize the overall signal intensities. Statistical significance within and between the groups was evaluated by Benjamini-Hochberg false discovery rate (FDR), as reported in the literature [30-32] FDR correction was implemented to control the spurious errors from multiple comparison tests and threshold of 0.05 and 0.2 was robust to call three prakriti specific differentially methylated regions. Quantile normalization and statistical test with FDR correction was performed using Gene Spring v12.6.

The differential methylated regions were identified at the probe level by both intra and inter-group analysis using relevant statistical measure. Intra-array probe level analysis was performed for the three prakriti using single sample t-statistics with Benjamini and Hochberg FDR correction to filter, most significant probes (p ≤ 0.05) with fold change of ≥1.5. Inter prakriti analysis was performed to find out the differentially methylated probes using analysis of variance with Tukey-post hoc test. Subsequently, differential methylated probes were identified with a fold change difference of ≥1.2 and 20% Benjamini and Hochberg FDR correction. The cross comparison of the intra-prakriti with inter-prakriti analysis was performed to find the prakriti specific differentially variable probes (mPSR). The differential methylated significant probes from the comparative analysis were hierarchically clustered using Pearson centered, Ward’s distance matrix. The identified mPSRs were analyzed across prakriti using Manhattan linkage and Ward’s distance matric. Intra-prakriti analysis by t-statistics with FDR correction and fold change analysis helped to find potential and, consistent methylated probes, whereas the ANOVA analysis identified differentially methylated probes. The identified mPSRs were analyzed with reference to CpG Islands and we designated differentially methylated if a) it contained at least one significant probe which fulfilled two criteria of probe level analysis and b) it should contain two or more methylated probes with fold change of ≥1.5 and p ≤ 0.05 in one of the three prakriti groups. For validation purpose, we have considered such differentially methylated CpG islands.

Association of differential mPSRs with genomic regions

The probe sequences were annotated using Agilent e-array technology file and Galaxy genome browser into different genomic regions based on the distance from TSS of the gene. The probes were mapped to 5′-UTR and 3′-UTR regions using RefSeq gene coordinates. The enrichment of mPSRs at chromosomes, CpG Island, CpG shores and different genomic position was statistically analyzed by Chi-square test (p < 0.05). To understand the association of histone modification and CpG methylation, we used EpiExplorer and the statistical significance was estimated by the Fischer Exact test with FDR correction [33]. Functional relevance for prakriti specific methylated gene sets was constructed by Gene Ontology (GO) term, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis by DAVID [34] and Funrich v2.1 (http://www.funrich.org/) tools with default parameters. The biologically enriched significant terms were deciphered by analyzing the prakriti specific methylated genes against the whole genome with hyper geometric test and p ≤ 0.05 without FDR correction was considered.

Validation by bisulfite sequencing (BS)

Bisulfite conversion of whole blood DNA was performed using the EZ DNA Methylation-Gold kit (Zymo Research, Orange, CA, USA) according to manufacturer’s protocol. In brief, about 1.5 μg of total blood genomic DNA for each dominant prakriti was treated sodium bisulfite and bisulfite specific PCR was performed. The amplicons were purified and sequenced using BigDye terminator cycle sequencing kit (ABI, USA) [35]. The quality of the sequencing results and absolute methylation at CpG position were analyzed using ESME software [36]. Two way ANOVA with Bonferroni multiple correction was used to account for the significance (p < 0.05). For each selected prakriti specific differentially methylated region, 20 dominant prakriti samples were chosen for validation. The primer list along with chromosomal coordinates and PCR conditions are provided in the supplementary information (Additional file 2: Table S1).

Results

Analysis of global methylation variation

Global DNA methylation estimation (5-mC) was performed for 36 samples randomly selected (n = 12 for each prakriti), using reverse phase HPLC methodology. We observed a significantly higher level in global methyl cytosine content (p < 0.05) in Vata when compared with Pitta and Kapha prakriti (Additional file 3: Figure S1).

Intra-prakriti methylation analysis

The intra-prakriti methylated probes were identified by providing a cutoff of fold change (FC) ≥1.5 with 5% FDR, based on the previous reports [30,37]. This analysis ensured the consistency of M-values across the prakriti groups and observed a total of 23,101, 23,276 and 21,014 significantly methylated probes in Kapha, Pitta and Vata prakriti respectively. Venn analysis showed a total of 2523, 3547 and 1535 methylated probes uniquely representing Kapha, Pitta and Vata prakriti respectively (Figure 1A). The distribution of probes based on their genomic location attributed to the nearest TSS site of a gene suggested a higher number of probes were associated with the gene body in Pitta prakriti as opposed to promoter associated probes in Kapha prakriti (Figure 1B). The Chi-square analysis on the distribution of methylated probes showed higher frequency at gene body, 5′UTR and 3′UTR, in case of Pitta (OR_p = 3.56, P < 0.001) as compared to other two prakritis. A significant decrease in promoter methylated CpG probes in all three prakriti (OR_k = 0.88, OR_p = 0.79 and OR_v = 0.86, p < 0.001) was also observed (Additional file 4: Figure S2).

Figure 1

Comprehensive differential DNA methylation analysis in prakriti. (A) Venn diagram for the identification of prakriti specific methylated probes. The numbers indicated in the black boxes represent the identified 501 mPSRs in which 23, 303 and 175 probes were found to be more methylated in Kapha, Pitta and Vata respectively. (B) Distribution of uniquely significant methylated probes with respect to CpG islands and CpG shores. The probes are classified on the basis of its position from TSS, as upstream promoter (−10Kb to -1Kb), promoter (−1Kb to transcription start site), 5′UTR, gene body (within the genic region), 3′UTR, Downstream (Up to 10Kb from the end of gene coordinates) and Intergenic. (C) Genome wide distribution of 501 mPSRs. Bar height indicate the degree to which each chromosome contains an unexpectedly high number of differentially methylated probes. Bar height is the –log10 values of fisher exact p-value derived from no FDR correction. The height of bars higher than black dashed line was considered to be significant (p ≤ 0.05). (D) Hierarchical cluster analysis of 501mPSR. M-values scale represents the relative methylation. Column represents the samples and rows represented by 501mPSRs. The arrays and probes belonging to Kapha, Pitta and Vata are shown in red, green and blue respectively.

We identified, 735 (639 genes) Pitta specific, 561 (509 genes) Kapha specific and 305 (262 genes) Vata specific methylated CpG islands, represented by two or more methylated probes respectively. Evaluation of multiprobe methylated CpG Islands in the prakriti across the genome with respect to extended promoter (−10000 bp to 500 bp from TSS) and gene body (inside the gene coordinates) showed variations in Kapha and Pitta prakriti (Additional file 5: Figure S3). In order to integrate fold change and P-value significance cut-off for the inclusion of multiprobe methylated CpG Island, a preliminary validation of frequently methylated CpG Island at the promoter of an imprinted gene, neuronatin (NNAT), was performed in all three prakriti samples (n = 20). Results showed an average methylation level of about 40-60% in individual CpG sites across all the samples from three different prakriti (Additional file 6: Figure S4).

It is reported that 70% of CpG content of whole genome is represented in CpG Islands and about 50% of CpG islands encompasses promoter regulatory element to facilitate the expression variations among the individuals [38]. We show, the distribution of LogR values to be less near TSS compared to upstream and downstream regions, suggesting that the majority of the CpG islands near promoter regions were unmethylated. A comparative analysis of methylated promoter associated genes with the whole blood gene expression analysis data (http://biogps.org/downloads/) [39] revealed an inverse correlation confirming our results (Additional file 7: Figure S5).

Differential methylation variation in prakriti

All 147 arrays irrespective of prakriti were subjected to quantile normalization. Analysis of variance (ANOVA) with Tukey post-hoc test at Benjamini Hochberg FDR of 20% correction was performed to enumerate differential methylated probes among prakriti samples. A total of 5447 differentially methylated probes were identified with fold change difference of ≥1.2 and FDR correction at P-value 0.2. These intra-prakrti probes were compared against significantly methylated probes of individual prakriti analysis. In all the, Kapha versus Pitta, Pitta versus Vata and Vata versus Kapha analysis, we observed 665 (Kapha vs Pitta), 755 (Pitta vs Vata) and 339 (Vata vs Kapha) probes to be differentially methylated. The supervised cluster analysis (Pearson’s distance metric and Ward’s linkage) of these differentially methylated probes (Additional file 8: Figure S6) indicated a distinct pattern of clusters of prakriti phenotypes. We determined a total of 501 prakriti specific differential methylated regions (mPSRs), represented by 303 probes overlapping with 269 CpG islands in Pitta prakriti, followed by Vata prakriti with 175 probes representing 139 CpG islands and finally Kapha prakriti with 23 probes and 19 CpG islands (Figure 1A) (Additional file 9: Table S2). The mPSRs were scattered across the genome and significantly (Chi-Square, p < 0.05) overrepresented at chromosome 3, 10 and 16 (Figure 1C). The supervised cluster analysis with Manhattan distance matrix and Ward’s linkage showed a prominent cluster of Pitta samples (Figure 1D) suggesting mPSRs could distinguish between prakriti. We did not find three prominent clusters of each prakriti due to a) selection of normal healthy individuals with > 60% of one prakriti and b) variation due to the influence of second dominant prakriti phenotype.

Further the mPSRs were significantly (p ≤ 0.05) enriched in the different genomic regions such as promoter and gene body associated with both CpG islands and shores (Additional file 10: Table S3). The crosstalk between DNA methylation and histone modifications illustrate influence of chromatin changes and has several implications including its role in normal variation, development, somatic cell reprogramming, and tumorigenesis [40]. To understand the possible variations, we compared the mPSRs with different regulatory regions using Epiexplorer tool [33]. We observed that the regions of H3K9ac, H3K27ac, H3K4me1, 2 and 3 and DNaseI hypersensitive regions were significantly (p≤0.05) augmented, underlining the effect of DNA methylation through intricate chromatin regulation. In addition to this, there was a significant (p ≤ 0.05) enrichment in the CTCF insulator sites (Additional file 11: Table S4).

Association of mPSRs with Lymphocyte and BMI

It is established that differences in cell types may also possess non-overlapping DNA methylation signatures within population [41]. Hence, we systematically tested for such regions by comparing the mPSRs with the additional lymphocyte microarray data. Lymphocyte DNA from 11 healthy individuals was analyzed in a similar way as that of intra-prakriti analysis. The unmethylated significant probes of lymphocyte were compared with methylatd significant probes of whole blood data and also with uniquely represented methylated probes. However, we found very few probes showing such variations among the prakriti. Further, it was interesting to note that the methylation values of 501 identified mPSRs of whole blood data were positively significantly correlated with lymphocyte data (Spearman correlation R = 0.15 and p = 0.0004), indicating that the identified prakriti specific regions were distinct and were not influenced by variation in the blood cell types tested (Additional file 12: Figure S7).

In a large population based study, we have previously shown an association between BMI and prakriti [23]. Hence, mPSRs were tested for their correlation to BMI by categorizing the hybridized samples into two groups: a) individuals with BMI less than 18 as low BMI (n = 23) and b) individuals with BMI greater than 25 as high BMI (n = 23). The arrays were quantile normalized and paired t-test with no correction method was used to identify differential methylated sequences. We found a total of 3989 probes with greater than ±1.2 fold change difference between high BMI and low BMI (p < 0.05) groups (Additional file 13; Figure S8). Among 501 mPSRs, 34 were common to differentially methylated BMI associated regions and were cluster distinctly with BMI as well as prakriti as well (Figure 2A and B). The genes corresponding to the 34 mPSRs were functionally associated with cellular adhesion or in the maintenance of cellular structural integrity and thus correlated with Kapha prakriti, characterized as energy of structure. Therefore cadherin22 precursor (CDH22) gene was selected for validation, to test the association of BMI group and prakriti.

Figure 2

DNA methylome analysis of prakriti and BMI (A) Comparative analysis of differential methylated probes of BMI and mPSRs (B) Hierarchical cluster analysis for the significant differential methylated probes commonly existing among prakriti and BMI. M-values were scaled such that red indicates higher methylation in a row and green indicates lower methylation. Clustering was performed using Ward’s hierarchical clustering algorithm with Pearson correlation distance as the distance metric. Rows correspond to samples and columns to probe regions.

Functional relevance

The functional relevance of the identified prakriti specific methylated genes was analyzed using DAVID [34]. The mPSR associated genes were analyzed and top significant enriched gene ontology and pathway terms with a p ≤ 0.2 were represented in table S5 (Additional file 14: Table S5). We anticipated that the observed differential methylated genes at different genomic regions may influence gene expression leading to phenotypic variation among prakriti. Although promoter methylation is known to show inverse correlation of expression, recent reports have suggested the significance of the gene body methylation and its role in gene expression regulation [42]. The rationale of the present study was to test the hypothesis for the prakriti specific methylation variation in the whole blood DNA. Considering these, we chose all the mPSRs genes and attempted to associate to different attributes as energy of structure, metabolism and motion for Kapha, Pitta and Vata respectively [11,43,44].

Several methylated sequences of Vata prakriti were found to be represented in diverse biological processes such as cell communication, transcription, various signal transduction pathways and embryo morphogenesis. The associated genes are involved especially in neuronal development (HOXB1, LMX1B, LHX1, LHX5, EN2; p ≤ 0.05). In addition, we observed association of NFIX gene with low BMI which is one of the characteristic features of Vata prakriti. Pitta prakriti showed enrichment of mPSRs for metabolism related Gene Ontology terms and pathways. The biological process such as regulation of hormone secretion (UNC13A; p = 0.01), regulation of nucleic acid metabolism (SKI, TP73, RUNX3, LMX1A, ZNF496, ZNF672, SOX11; p = 0.03) and several others were significantly enriched. Apart from these, signal transduction pathways such as electron transport system, NOTCH, SIP, AKT, BCR signaling associated genes were also enriched in Pitta prakriti. Interestingly, we have observed regulation of OCA2, MC1R genes, which are associated with white to blonde or red hair colour in Pitta phenotype. Kapha mPSR associated genes were significantly enriched in cell growth/maintenance (PRELP, COL2A1, ACTR10, CDH22; p = 0.03), cytoskeleton anchoring activity (PRELP; p = 0.03) and cellular adhesion (CDH22, CADM1; p = 0.05) function.

Microarray validation

Experimental validation on prakriti specific hypermethylated genes such as LHX1, SOX11 and CDH22 was performed. LHX1- a LIM homeobox 1, is a member of LIM domain which has cysteine-rich zinc-binding domain function as transcription regulator. The gene essentially involved in the development of retinal, urogenital system and localization and differentiation of tissues [45]. Recent report highlighted the disorganized circadian rhythm due to LHX1 expression deficiency in suprachiasmatic nucleus [46]. Individuals with Vata prakriti have irregular sleeping and eating habits which is in accordance with the fragile circadian rhythm, and may be associated by LHX1 promoter methylation. The average methylation of the amplicon was found to be more in Vata in comparison to Pitta and Kapha. The post-hoc analysis with Bonferroni correction revealed a significant (p < 0.05) higher methylation at CpG site 47, 55 and 65 in Vata prakriti (Additional file 15: Figure S9).

Pitta methylated, SOX11, an intron less gene encodes a member of the SOX (SRY-related HMG box) family of transcription factors mainly involved in the regulation of embryonic development and in the determination of the cell fate, present in the chromosome 2p25.2. The protein encoded is a transcriptional regulator possessing diverse functions namely, visual [47], renal system development and in tumorigenesis [48]. Individuals with Pitta prakriti governs body metabolism, temperature regulation and development/proper function of eyes. Although, these provide only a causal relationship between the methylation and phenotype, a functional experimental study is essential. We observed that overall average methylation at all CpG sites was found to be higher in the Pitta prakriti when compared to other two (Additional file 16: Figure S10).

A Kapha specific and BMI associated cadherin family protein, cadherin 22 precursors (CDH22) was validated. The gene is mainly involved in the morphogenesis and tissue formation in neural and non-neural developmental process. The recent reports suggested that several known SNPs in the CDH22 were associated with type II diabetes and higher BMI [49]. A differentially methylated probe’s M-value is significantly positively correlated to the BMI (Spearmann correlation R = 0.22 and p < 0.05) (Figure 3A). The CDH22 amplicon was 301 bp in length, present at the 5′UTR of the gene, encompassing a total of 40 CpG sites out of which 33 CpG sites were analyzed in all the three prakriti. The fold change difference between prakriti underlined a higher methylation level for Kapha prakriti which was in concordant with the direct BS (Figure 3B and C). Overall, amplicon methylation was significantly (p < 0.05) high in Kapha and high BMI when compared to other prakriti and low BMI respectively. Apparently, absolute methylation of at individual CpG sites of individuals with Kapha prakriti and high BMI was found to be similar as opposed to individuals with Vata prakriti who have low BMI (Figure 3D and E).

Figure 3

DNA methylation profiling at cadherin 22, precursor (CDH22) in prakriti and BMI. (A) Scatter plot illustrating positive, significant correlation (R = 0.22 and p ≤ 0.05) between the methylation values of CDH22 and BMI. (B) Distribution of log fold change differences of CDH22, 5′UTR associated CpG island. The differences with positive values in yellow and negative with blue indicate level of methylation variation between prakriti. The height of the individual bars represents the extent of methylation differences. Shaded in blue is a probe region showing a significant difference at p ≤ 0.2 between prakriti. The offset scale for log2 fold differences was set to −0.2 to 0.2. (C) Bisulfite genome sequence analysis. Methylation cut off >15% was considered to call methylated sites and are represented in black circles. The circles in white were considered as unmethylated with respect to the cut-off chosen. Rows indicate the sample and columns represent individual CpG site. (D) and (E), is the cumulative CpG site methylation variation in prakriti and BMI respectively. Each circle represents the overall methylation of amplicon for an individual. The methylation values in the groups is indicated with mean ± standard error and *refer to the significance (*for p ≤ 0.05 and **for p ≤ 0.01).

Discussion

DNA methylation plays central role in gene regulation and may influence variations in human phenotypes. Although DNA methylation is reported to be associated with several complex diseases including cancer [50,51], only limited studies have demonstrated its association with variations in the normal phenotypes. Ayurveda, an ancient Indian traditional medicine, classifies human population into three distinct prakriti phenotype as Vata, Pitta and Kapha prakriti [10,13]. The prakriti is a relative proportion of enumerated characteristics and may remain constant throughout the life span of an individual [52]. Ethnicity, familial characteristics, and habitat appear to influence phenotypic variability through their effect on prakriti. Similar factors also appear to contribute to inter-individual variability at the genetic or epigenetic levels [53].

However, precise molecular basis for the manifestation of such an effect remains to be understood. In the present study, we made an attempt to identify the methylation signatures associated with different prakritis and the potential impact on determining distinct prakriti specific phenotypes. Our study showed that large proportions of DNA methylation patterns are common between prakritis; observed differences in methylation signatures between prakriti suggests the presence of different mechanisms that may influence and sustain the expression of key regulator genes involved in the manifestation of distinct phenotypes.

From, intra-prakriti analysis, Pitta prakriti showed higher number of methylated CpG islands in gene body region, while promoters regions were more methylated in Kapha prakriti. CpG island associated within and outside the promoters as well as in the gene body can directly or indirectly affect the rate of transcription [3,54]. Methylation profile within individual CpG island is random and stochastic [55], therefore averaging the M-values or Z-score values for each CpG islands may not be sufficient to identify the differentially methylated CpG islands. Hence we used a simple methodology to identify differentially methylated signatures in the prakriti by considering methylated multiple probes associated CpG islands represented by a defined number of differentially methylated probes. Further, fold change cut-off of ≥1.5 for the methylated probe was used and as a proof of concept, a commonly methylated promoter region of imprinted NNAT gene was validated. To identify differential methylation between prakriti, ≥1.2 fold change cut-off, with FDR correction of 0.2 was provided as reported earlier [30]. We have also validated promoters and 5′UTR of three genes-LHX1, SOX11 and CDH22 which are found to be hypermethylated in Vata, Pitta and Kapha prakriti respectively. Validation of LHX1, SOX11 and CDH22 in a panel of prakriti samples showed concordance with the findings from microarray data.

Differences in DNA methylation patterns have been reported among various cell types in the human peripheral blood [41]. The comparative analysis of DNA methylation profile of lymphocyte and total blood DNA in specific to prakriti were minimal and do not influence much on the mPSRs. Thus, in this preliminary study, the 501 different mPSRs identified were sufficient to classify the individuals based on the prakriti; thereby providing a molecular evidence for the ancient concept of Ayurveda. The distribution of differentially methylated regions between the prakriti as well as within the prakriti were not uniform, suggesting that stochastic DNA methylation differences in prakriti may not be limited for the short-listed genes, but may influence multiple functional gene networks in multiple chromosomes. Recent reports have suggested that gene body DNA methylation may be involved in the regulation of gene expression either directly or indirectly [3]. The function of gene body DNA methylation remains unclear and recent studies have reported the potential for how gene body DNA methylation affects gene expression [56,57]. The gene-body-specific methylation has also been reported to be associated to replication timing and histone modification [58,59]. Our results have shown the enrichment of methylated CpG islands in gene body regions for Pitta and promoter regions for Kapha prakriti. This suggests that complex epigenomic variations within prakriti are regulated by genomic position of CpG methylated sequences.

Unlike prakriti, DNA methylation reprogramming events may also be determined during gestation, development and pregnancy. Several reports suggested that DNA methylation of cord, placental and whole blood were influenced by the season of conception, maternal diet, starvation, maternal stress by depression and anxiety, smoking and socio-economic conditions. Some of these have been implicated in stable phenotypes like low birth weight, height and obesity [60,61]. Therefore the identified mPSRs were correlated with BMI and showed the association of Kapha prakriti with high BMI and Vata with low BMI, which is in consistency with our earlier report [23]. Global 5-mC content between Kapha and Pitta there was no significant difference indicating that such similarity may be limited at the global level. Our study shows the importance of both global and gene specific methylation associated with prakriti phenotypes.

Gene ontology and pathway analysis for the mPSRs suggests the enrichment of gene related to characteristics of specific prakriti. The Kapha specific methylated genes were enriched in cellular adherence, protein localization and transport. The Pitta associated methylated specific genes were enriched in the various transcriptional DNA binding motifs, developmental morphogenesis and B-cell, T-cell receptor signaling pathway. Vata specific methylated genes were found to be associated with DNA sequence specific binding proteins involved in transcription regulation, nervous system development, metal ion binding motifs regulating cellular motion and cell-cell signaling. Thus the identified methylated genes may either directly or indirectly have an influence on the prakriti phenotypes. For example, Kapha specific mPSRs genes found enriched in cell adhesion molecules (CDH22, CADM1, CDH11) (p<0.05), there by regulate structure and integrity of the cells. Further the concordance of CDH22 methylation validation was evident to find the association of Prakriti and BMI. The mechanisms and causes which sustain these changes in relation to prakriti determination needs further investigation and also demands for the case–control or retrospective large scale epigenome wide studies based on prakriti. However, we could able to demonstrate presence of the methylation signatures correlating to prakriti and its distinct phenotype specifically the BMI as an example. Thus, our results provide evidence for the identified genes for the specific prakriti attribute directly or indirectly to bring about defined phenotypes and suggest the impact of DNA methylation on specific prakriti.

Conclusion

Therapy in Ayurveda is based on the premise of maintaining human body or any living cell as an open system, which can partake with the external environment. Prakriti may govern functions at cellular, tissue and at system levels maintaining equilibrium determined by genes which in turn are regulated by epigenetic cross talk especially through DNA methylation. These may provide basis to explore whether the DNA methylation programmed at the time of birth of an individual that will mimic the prakriti constitutions. This is suggested on the basis of genetic, biochemical, hematological or anatomical features already identified for the specific prakriti. Therefore, our study on the epigenetic basis of prakriti is a step towards deciphering diagnostic and prognostic utilities of traditional medicine by identifying genome-wide DNA methylation variants. In summary, we have identified DNA methylation signatures that distinguish the three major prakriti, which are responsible for imparting variations in humans. Taken together, DNA methylation is probably coupled to chromatin regulation as a contributor to different prakriti phenotypes such as Vata, Pitta and Kapha. The differentially methylated regions identified in our study needs further investigation as they may provide insight into epigenotype–phenotype correlation and that underlining molecular mechanisms of prakriti manifestations.