Epigenetics is a study on genetic modifications, which occur in a cell without any change in the DNA sequence. Epigenetic regulation influences gene expression by turning its switch on or off. They contribute mainly to DNA methylation, histone modifications and RNA processing. The role of epigenetic regulation is widely studied in various disciplines. Its role in abiotic stress tolerance is a challenge for plant researchers.
Drought is one of the major abiotic stresses and an important limiting factor of rain-fed agriculture. Evidences on genome wide methylation patterns suggest hypo  or hyper-methylation  as response to stress and thus regulating gene expression. Recently, with the advent of next generation sequencing, DNA methylation is detected using sodium bisulfite conversion of cytosines to uracils, where 5’-methylcytosines remain unaltered. Bisulfite-treated DNA is then, amplified and sequenced . Genome wide methylation studies on various model plants are headway to address drought tolerance but there is no emphasis on the organellar genomes.
Chloroplast and mitochondria are crucial support to energy metabolism occurring in a plant. They respond in an effective, energy-efficient way to combat continuous changes in the environment. They play a significant role in cellular metabolism and therefore, are expected to be first in stress recognition . In the year 1998, Simkova  had studied the methylation of Mitochondrial DNA in carrot using Southern hybridization and in 2012. Huang et al.  reported their study on mitochondrial and chloroplast DNA in the woody plant, Sequoia sempervirens using HPLC. In the recent years, tremendous data have been generated with respect to mitochondrial DNA methylation in humans. Mawlood et al. (2016)  suggest mtDNA methylation as biological marker for forensic age-prediction and health status measurement. A Study by Sun et al., (2018)  on tumorigenesis highlight the influences that the nuclear and mitochondrial genomes have in setting mtDNA methylation patterns to regulate mtDNA copy number. Another study on distribution and dynamics of MtDNA methylation suggests its role during gametogenesis . Similarly, methylation patterns in mitochondrial DNA are related to cardiovascular disease , Down’s syndrome , obesity  etc. However, limited work has been done on DNA methylation in plants as compared to Humans in general and mtDNA methylation in particular. Feng et al. (2016)  have studied cadmium exposed DNA methylation pattern in rice. Another study by Kawakatsu et al. (2017)  involves comparison of time-series methylomes of dry and germinating seeds to publicly available seed development methylomes and has concluded the role of DNA methylation in seed dormancy. Role of DNA methylation is seen in fruit ripening  and abiotic stresses .
However, mitochondrial and chloroplast DNA methylation has been faintly studied in plants. Very recently, Muniandy et al., (2019)  have published their comparative study on methylation of chloroplast and amyloplast genomes of rice.
Therefore, in this study we have worked on the DNA methylation patterns in mitochondria and chloroplast of maize. We have been working on transcriptomes of Zea mays Z59 for drought tolerance in crop growing seasons, Kharif and Rabi  . In this study, we have worked on the methylation pattern of mitochondria and chloroplast and their relation with the gene expression (transcriptome data of the same cultivar) in response to abiotic stress, drought.
2. Materials & Methods
2.1. Plant Material & Growth Conditions
The seeds of Zea mays Z59 were sown in pots with uniform soil under greenhouse conditions. The growth conditions, viz., temperature, rainfall, relative humidity and other parameters were recorded for both the dry and wet seasons (Supplementary File 2). Pots were well-watered until the reproductive stage and drought was imposed from this stage till maturity by withholding irrigation for half of the plants (16 out of 32). After the onset of drought, the leaves were dark adapted for 15 mins and the readings were measured as OJIP parameters and Performance Index, PI(abs) using portable fluorometer, Handy PEA+ (Hansatech instruments ltd., Norfolk) at regular intervals. Performance Index (PI) includes three independent parameters: 1) density of fully active reaction centers (RCs); 2) efficiency of electron movement by trapped exciton into the electron transport chain beyond the QA; and 3) the probability that an absorbed photon will be trapped by RCs. PI(abs) is calculated by the formula
where: F0 means fluorescence intensity at 50 μs, FJ is fluorescence intensity at the J step (at 2 ms), FM represents maximal fluorescence intensity, VJ is relative variable fluorescence at 2 ms calculated as , M0 represents initial slope of fluorescence kinetics, which can be derived from the equation:
The extent of drought stress was observed by measuring relative water content and soil volumetric water content. The relative water content was measured from the third leaf at regular intervals and was calculated by the formula  (fresh weight − dry weight)/(turgid weight − dry weight) × 100. The soil volumetric water content was measured with Field Scout TDR 350 probe (Specturm Technologies Inc.). DNA was isolated from the leaf tissue (lamina, topmost leaf) from both the sets of irrigated and water-stressed plants and was subjected to bisulfite sequencing to study the methylation pattern during water stress.
2.2. Library Preparation and Bisulfite Sequencing
Libraries were prepared using 200 ng of fragmented DNA (~200 bp) using NEBNext Ultra DNA Library Prep Kit for Illumina (NEB), according to the manufacturer’s protocol. DNA is then, treated with sodium bisulfite followed by PCR amplification. Unmethylatedcytosines are deaminated to uracil, and read as thymine after PCR amplification and sequencing, whereas methylated cytosines are not converted and therefore, are read as cytosine after PCR amplification and sequencing  .
2.3. Assembly of the Bisulfite Treated Sequencing Reads
The raw reads were quality checked using FASTQC  and the adapters were trimmed using cutadapt . The clean reads were aligned to the B73 reference genome (AGPv3) using the Bismark aligner  using default parameters. Methylated cytosines were extracted from aligned reads using the Bismark methylation extractor available in the same package. The percent methylation in each CpG was calculated by the formula (number of reads with methylated C/total reads) × 100.
2.4. Sliding Window Approach
The data was sorted chromosome-wise and the methylation data of the organellar genomes were used in further downstream analysis. A sliding window of window size 1000 bp and sliding size 100 bp were used to generate methylation count for each window. The peaks showed the hyper-methylated and hypo-methylated regions.
2.5. Methylation Data Analysis
The methylation data in the CpG context of the organellar chromosomes were taken from the total genome wide methylation data. The regions of methylation were manually mapped to the genome feature file of the reference genome for both Mitochondria and Chloroplast for annotation of the methylated regions. Gene-wise methylation count was thus generated giving information of the genes that are hyper and hypo methylated due to stress.
2.6. Co-Relation with RNA-Seq Data
RNA-seq reads are available for Zea mays Z59 . The clean RNA-seq reads were aligned to the organellar chromosomes of B73 reference genome (AGPv3) using Bowtie2 aligner . The mapped reads for both Mitochondrion and Chloroplast were de novo assembled separately using Trinity  with K-mer = 25 for both irrigated and water stressed plants. The abundance of the trinity generated transcripts was calculated using RSEM . Normalization of the expression values of irrigated and water stressed transcripts of both mitochondria and chloroplast was generated separately as TPM values considering the length of the transcript and the number of reads that mapped to the transcript. The transcripts were annotated using blastx  and GeSeq . Circular maps for both mitochondria and Chloroplast Zea mays Z59 were generated using OGDRAW  after ordered contigs were obtained from Progressive Mauve alignment . For each gene, the expression values were plotted against the methylation count to study the co-relation of gene expression and methylation. Statistical analysis of the expression and methylation data was performed using SAS .
3.1. Drought Stress and Photosynthesis
The plants were monitored at regular intervals. The sliding down of relative water content and soil water volumetric content readings shows the onset of drought in water-stressed plants. The OJIP parameters and performance index, PI(abs) did not show significant difference between irrigated and water-stressed samples initially, later which there was a slight decline in water-stressed plants compared to irrigated as the drought period prolonged (Figure 1). This indicated that photosynthesis may be affected after a prolonged period of drought in this drought tolerant cultivar.
3.2. Bisulfite Illumina Sequencing and Analysis
Whole genome bisulfite sequencing of Zea mays Z59 resulted in 260485668 and 204847228 number of paired-end reads in irrigated and water-stressed samples. Assembly of the reads with the reference genome B73 (AGPv3) with Bismark
Figure 1. (a) Chlorophyll a fluorescence transients in irrigated and water stressed samples at 10 (Blue, Violet), 30 (Grey, Black), 50 (Pink, Green) days of drought respectively; (b) PI(abs) values at regular intervals of drought; (c) & (d) Soil volumetric and relative water content in irrigated and water-stressed samples at regular intervals of drought.
Aligner resulted in the identification of genome-wide methylated CpG sites. The methylated base pair position along with its tri nucleotide context is generated. The methylation data of the organellar genomes were used for further analysis (Supplementary File 1). From, the overall methylation count, water-stressed sample shows hypo methylation compared to irrigated in both Mitochondria and Chloroplast. Classification of methylated positions based on gene type shows highest methylation count in protein-coding genes in both Mt and Pt. However, in Mt, good number of rRNA, tRNA and pseudogenes also show methylation as compared to chloroplast where anti-sense RNA showed methylation (Figure 2).
With window size 1000 bp and sliding size 100 bp, the methylation count of the organellar genomes showed highest peak between 419,800 to 420,800 base pair region in Mitochondria (Figure 3) and 36,900 to 37,900 base pair region in chloroplast in both irrigated and water-stressed samples (Figure 4). Annotation of the methylated regions showed high methylation in nad4 gene in Mitochondria and psaB, atpA, psaA, psbD, psbB, ropB, rpoC1 genes in Chloroplast.
3.3. Co-Relation of Methylation Data with RNA-Seq Data
The RNA-seq reads of Zea mays Z59 were aligned to the organellar B73 reference genome using bowtie 2 . De novo assembly of the mapped reads generated 259 transcripts in mitochondria having GC content 44.94 and 59 in chloroplast having GC content 36.68. The normalized expression values of these transcripts show high expression of rps3, rps2A, ccmFC, atp1 and many uncharacterized genes in mitochondria and psbA, psbD, psbc, psaA, atpA, ndhk, ndhB1 & B2 genes in chloroplast. The methylation count of these genes were plotted against the expression values (Figure 5) and it is typically observed that the genes that are hypomethylated in water stressed plants when compared to irrigated plants are up-regulated, and down-regulated genes are hypermethylated. This was also confirmed by the statistical analysis where a negative correlation was observed (Figure 6). The genes are related to energy metabolism in mitochondria and photosystem I & II genes in Chloroplast. Here, we understand that
Figure 2. Methylation count in different gene types.
Figure 3. Sliding window graph of methylated sites in mitochondria showing row number of the window size on the X-axis and methylation count on the Y-axis. X-axis Unit is basepair (bp) Y-axis Unit is Number of methylated sites.
Figure 4. Sliding window graph of methylated sites in chloroplast showing row number of the window size on the X-axis and methylation count on the Y-axis. X-axis Unit is basepair (bp) Y-axis Unit is Number of methylated sites.
Figure 5. Gene expression vs methylation count in mitochondria and chloroplast. X-Axis Gene accession number, Y-Axis TPM, Z-Axis Number of methylated sites.
Figure 6. Matrix plot of Chr Mt and Pt showing negative correlation between gene expression and methylation count. X-Axis: TPM, Y-Axis: Number of methylated sites.
due to hypomethylation of these stress responsive genes, they are up regulated indicating maintenance of photosystem II performance even in water stressed conditions. The annotated circular plots of Zea mays Z59 were built for both chloroplast and mitochondria (Figure 7 & Figure 8).
Figure 7. Circular plot of chloroplast genome of Zea mays Z59.
3.4. Data Availability
The data is submitted to NCBI and can be accessed using the following accession numbers.
Figure 8. Circular plot of mitochondrial genome of Zea mays Z59.
Plant vitality during water stress is an integrative study, involving network of genes and their relationships. The organellar genomes of the plants also play an essential role like the nuclear genome during abiotic stress. Our study on methylation and gene expression patterns of the organellar genes during water stress in Maize has given wealth of information for plant breeders about this drought tolerant genotype Z59.
The relationship between methylation, gene expression and PI(abs) & OJIP parameters is like a triangular hypothesis. DNA methylation affects gene expression which in turn monitors PI(abs) & OJIP readings. Hypo methylation of genes in water stressed samples showed up regulation of the same genes, thereby showing tolerance to water stress which is reflected in the PI(abs) and the OJIP readings.
When leaves are dark-adapted for a minimum of 15 minutes all photochemistry in the chloroplasts stops, on subsequent exposure to light, photochemistry starts and chlorophyll fluorescence rises from a low to a high. This fast fluorescence transient rise is directly related to the Photosystem II activity. OJIP transient rise reflect the reduction reactions taking place in the PSII which involves QA, one electron accepting bound PQ, QB and two electron accepting mobile PQ. When all the reaction centres are open upon light interception it is termed as level O where there is no reduction of QA, here the fluorescence intensity lasts for 10ms. After this the fluorescence transient rises from O to J which is when there is the net reduction of QA which is the stable primary electron acceptor of PS II to Q− A, this rise lasts for 2 ms. The J to I phase is due to the resultant reduced state of closed Reaction Centres (RCs) QA− QB−, QA QB2− and QA− QB H2 which lasts for 2 - 30 ms. At 300 ms the plastoquinol pool is maximally reduced and this along with maximum concentration of QA−QB2 is represented by the level P in the OJIP curve .
The quantitative information on the plant performance under stress was studied using PI(abs). It reveals the functionality of both photosystem I & II . Both PI(abs) and OJIP transient values are observed akin in water stressed samples as in irrigated samples probably because of the upregulation of photosystem II genes in water stressed samples. Up regulation of functional genes during stress can be attributed as one of the mechanisms to combat stress in drought tolerant genotypes like Zea mays Z59.
Joel (2013)  speculated hyper-methylation as an indicator of drought susceptibility and hypo-methylation for drought tolerance. We have also observed hypomethylation, overall, in water stressed sample for both Mitochondria and Chloroplast indicating up regulation of genes like rps3, rps2A, ccmFC, atp1 and many uncharacterized genes in mitochondria and psbA, psbD, psbC, psaA, atpA, ndhk, ndhB1 & B2 genes in chloroplast. Ribosomal protein S3 and S2 (rpS3 & rps2) are known to play critical roles in ribosome biogenesis and DNA repair during increase of cellular ROS levels as response to stress . Cytochrome c biogenesis (ccmFC) gene is known to form a complex with CCMH, CCMFN1 and CCMFN2 that performs the assembly of heme with c-type apocytochromes in mitochondria. According to the study by Evans et al., (2019) , thisgene is essential for mitochondrial function. Disruption of this might lead the mitochondrion to be non-functional. Mitochondrial ATP1 and atpA of chloroplast produces ATP from ADP and ATP is an immediately available energy form. The up regulation of this gene indicates enhancement of energy requirements as response to stress .
In chloroplast, the upregulated genes are psbA, psbD, psbC of photosystem II and psaA of photosystem I. The genes ndhk, ndhB1 & B2 are subunits of the chloroplast NAD (P) H dehydrogenase (NDH) complex involved in photosystem I cyclic electron transport and chlororespiration . This indicates that the photosystem I and II activity remains unaltered during water stress.
Whole genome bisulfite sequencing of irrigated and water stressed Zea mays Z59 plants and assembly of the reads on to the reference organellar genomes showed hypomethylation in water stressed plants when compared to irrigated plants. Annotation of the methylated genomes showed that genes related to photosystem I & II in chloroplast and nad 4 gene in mitochondria were hypo methylated in the water stressed sample; RNA-seq analysis of transcriptomic reads mapped to the same reference showed regulation of hypomethylated genes in water stressed plants.
Ms. B. Divya Bhanu acknowledges Department of Science and Technology, Government of India for financial support vide reference no SR/WOS-A/LS-597/2013 under Women Scientist (A) scheme to carry out this work. We acknowledge ICAR-CRIDA and National Innovations in Climate Resilient Agriculture (NICRA) for the facility provided for growth of maize plants. We also acknowledge CIMMYT-Hyderabad regional centre for sharing with us the seeds of inbred line Z59 with drought tolerance characters.
Mt—Mitochondria, Pt—Chloroplast, PI—Performance Index, Abs—Absolute, RC—Reaction Centres, PCR—Polymerase Chain Reaction, NCBI—National Centre for Biotechnology Information.
Supplementary File 1
Methylation results of chloroplast and mitochondria of irrigated and water-stressed samples https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/EMIFGY.
Supplementary File 2
 Wada, Y., Miyamoto, K., Kusano, T. and Sano, H. (2004) Association between Up-Regulation of Stress-Responsive Genes and Hypomethylation of Genomic DNA in Tobacco Plants. Molecular Genetics and Genomics, 271, 658-666.
 Rico, L., Ogaya, R., Barbeta, A. and Penuelas, J. (2014) Changes in DNA Methylation Fingerprint of Quercus ilex Trees in Response to Experimental Field Drought Simulating Projected Climate Change. Plant Biology, 16, 419-427.
 Van, A.O., Clercq, I.D., Ivanova, A., Law, S.R., Van, B.F., Millar, A.H. and Whelan, J. (2016) Mitochondrial and Chloroplast Stress Responses Are Modulated in Distinct Touch and Chemical Inhibition Phases. Plant Physiology, 171, 2150-2165.
 Huang, L.C., et al. (2012) DNA Methylation and Genome Rearrangement Characteristics of Phase Change in Cultured Shoots of Sequoia sempervirens. Physiologia Plantarum, 145, 360-368.
 Mawlood, S.K. Dennany, L., Watson, N., Dempster, J. and Pickard, B.S. (2016) Quantification of Global Mitochondrial DNA Methylation Levels and Inverse Correlation with Age at Two CpG Sites. Aging, 8, 636-641.
 Sun, X., Vaghjiani, V., Jayasekara, W.S.N., Cain, J.E. and John, J.C.S. (2018) The Degree of Mitochondrial DNA Methylation in Tumor Models of Glioblastoma and Osteosarcoma. Clinical Epigenetics, 10, Article No. 157.
 Sirard, M.A. (2019) Distribution and Dynamics of Mitochondrial DNA Methylation in Oocytes, Embryos and Granulosa Cells. Scientific Reports, 9, Article No. 11937.
 Baccarelli, A.A. and Byun, H.M. (2015) Platelet Mitochondrial DNA Methylation: A Potential New Marker of Cardiovascular Disease. Clinical Epigenetics, 7, Article No. 44.
 Infantino, V., et al. (2011) Impairment of Methyl Cycle Affects Mitochondrial Methyl Availability and Glutathione Level in Down’s Syndrome. Molecular Genetics and Metabolism, 102, 378-382.
 Feng, S.J., et al. (2016) Variation of DNA Methylation Patterns Associated with Gene Expression in Rice (Oryza sativa) Exposed to Cadmium. Plant, Cell & Environment, 39, 2629-2649.
 Liu, R.A. (2017) DEMETER-Like DNA Demethylase Governs Tomato Fruit Ripening. Proceedings of the National Academy of Sciences of the United States of America, 112, 10804-10809.
 Muniandy, K., Tan, M.H., Song, B.K., Ayub, Q. and Rahman, S. (2019) Comparative Sequence and Methylation Analysis of Chloroplast and Amyloplast Genomes from Rice. Plant Molecular Biology, 100, 33-46.
 Divya Bhanu, B., Ulaganathan, K., Shanker, A.K. and Desai, S. (2016) RNA-Seq Analysis of Irrigated vs. Water Stressed Transcriptomes of Zea mays Cultivar Z59. Frontiers in Plant Science, 7, 239.
 Bhanu, B.D., Ulaganathan, K. and Shanker, A.K. (2019) Seasonal Variation in Expression Pattern of Genes in Irrigated and Water Stressed Transcriptomes of Zea mays Z59. Journal of Plant Biochemistry and Biotechnology, 28, 271-279.
 Pieczynski, M., et al. (2013) Down-Regulation of CBP 80 Gene Expression as a Strategy to Engineer a Drought-Tolerant Potato. Plant Biotechnology Journal, 11, 459-469.
 Li, Y. and Tollefsbol, T.O. (2011) DNA Methylation Detection: Bisulfite Genomic Sequencing Analysis. In: Tollefsbol, T., Ed., Epigenetics Protocols, Humana Press, Totowa, NJ, 11-21.
 Landin, D.B.V., Pflüger, J. and Lister, R. (2018) Generation of Whole Genome Bisulfite Sequencing Libraries for Comprehensive DNA Methylome Analysis. In: Jeltsch, A. and Rots, M., Eds., Epigenome Editing, Humana Press, Totowa, NJ, 291-298.
 Krueger, F. and Andrews, S.R. (2011) Bismark: A Flexible Aligner and Methylation Caller for Bisulfite-Seq Applications. Bioinformatics, 27, 1571-1572.
 Li, B. and Dewey, C.N. (2011) RSEM: Accurate Transcript Quantification from RNA-Seq Data with or without a Reference Genome. BMC Bioinformatics, 12, Article No. 323.
 Conesa, A., et al. (2005) Blast2GO: A Universal Tool for Annotation, Visualization and Analysis in Functional Genomics Research. Bioinformatics, 21, 3674-3676.
 Greiner, S., Lehwark, P. and Bock, R. (2019) Organellar Genome DRAW (OGDRAW) Version 1.3.1: Expanded Toolkit for the Graphical Visualization of Organellar Genomes. BioRxiv, 545509.
 Darling, A.E., Mau, B. and Perna, N.T. (2010) ProgressiveMauve: Multiple Genome Alignment with Gene Gain, Loss and Rearrangement. PLoS ONE, 6, e11147.
 Stokes, M., Chen, F. and Gunes, F. (2014) An Introduction to Bayesian Analysis with SAS/STAT® Software. Proceedings of the SAS Global Forum 2014 Conference, SAS Institute Inc, Cary, USA.
 Zhu, X.G., et al. (2005) Chlorophyll a Fluorescence Induction Kinetics in Leaves Predicted from a Model Describing Each Discrete Step of Excitation Energy and Electron Transfer Associated with Photosystem II. Planta, 223, 114-133.
 Strasser, R.J., Tsimilli-Michael, M. and Srivastava, A. (2004) Analysis of the Chlorophyll a Fluorescence Transient. In: Papageorgiou, G.C. and Govindjee, Eds., Chlorophyll a Fluorescence. Advances in Photosynthesis and Respiration, Vol. 19, Springer, Dordrecht, 321-362.
 Gayacharan and Joel, A.J. (2013) Epigenetic Responses to Drought Stress in Rice (Oryza sativa L.). Physiology and Molecular Biology of Plants, 19, 379-387.
 Kim, Y., Kim, H.D. and Kim, J. (2013) Cytoplasmic Ribosomal Protein S3 (rpS3) Plays a Pivotal Role in Mitochondrial DNA Damage Surveillance. Biochimica et Biophysica Acta (BBA), 1833, 2943-2952.
 Evans, D.L., Hlongwane, T.T., Joshi, S.V. and Pachón, D.M.R. (2019) The Sugarcane Mitochondrial Genome: Assembly, Phylogenetics and Transcriptomics. PeerJ, 7, e7558.
 Kosová, K., Vítámvás, P., Urban, M.O., Prásil, I.T. and Renaut, J. (2018) Plant Abiotic Stress Proteomics: The Major Factors Determining Alterations in Cellular Proteome. Frontiers in Plant Science, 9, 122.
 Peng, L., Yamamoto, H. and Shikanai, T. (2011) Structure and Biogenesis of the Chloroplast NAD(P)H Dehydrogenase Complex. Biochimica et Biophysica Acta (BBA), 1807, 945-953.