Identification, methylation profiling, and expression analysis of stress-responsive cytochrome P450 genes in rice under abiotic and phytohormones stresses

ABSTRACT The cytochrome P450 (CYP) is a large and complex eukaryotic gene superfamily with enzymatic activities involved in several physiological and regulatory processes. As an objective, an in-silico genome-wide DNA methylation (5mC) analysis was performed in rice (Oryza sativa cv. Zhonghua11), and the epigenetic role of CYPs in two abiotic stresses was observed. Being a stable representative mark, DNA-methylation alters the gene expression under stressful environmental conditions. Rice plants under salinity and drought stresses were analyzed through MeDIP-chip hybridization, and 14 unique genes of the CYP family were identified in the rice genome with varying degrees of methylation. The gene structure, promoter sequences, and phylogenetic analysis were performed. Furthermore, the responses of CYPs to various abiotic stresses, including salinity, drought, and cold were revealed. Similarly, the expression profile of potential CYPs was also investigated under various phytohormone stresses, which revealed the potential involvement of CYPs to hormone regulations. Overall, the current study provides evidence for CYP’s stress regulation and fundamental for further characterization and understanding their epigenetic roles in gene expression regulation and environmental stress regulation in higher plants.


Introduction
Global climate changes exert a plethora of abiotic stresses on plants, including salinity and drought that ultimately lower crop production. 1,2 Plants have evolved a variety of physiological and biochemical mechanisms to cope with environmental hazards. However, genomics and molecular biology approaches revealed numerous stress-responsive genes to tolerate the related hazards. 3 Meanwhile, altering the stress-responsive genes expression level with their protein level helps to understand the plant's ability to tackle adverse environmental conditions. The Cytochrome P450 (CYP) is an ancient and ubiquitous enzyme protein that is mainly accountable for the oxidative metabolism of diverse exogenous and endogenous compounds. 4 The CYPs enzymes are designated due to their ability to absorb the light at 450 nm and catalyzed numerous xenobiotic and biosynthetic pathways such as xenobiotics detoxification, carbon assimilation, secondary metabolites biosynthesis, drug deactivation, defense-related compound, UV protection, phytohormones, pigmentation, and procarcinogen activation. [5][6][7][8][9] Nomenclature and classification of CYPs are based on peptide similarity. According to this system, two CYP proteins sharing a sequence identity of greater than 40% belong to the A-type family, and 55% belong to the non-A-type subfamily. In plants, the A-type was included in CYP71 clads, whereas the non-A type was grouped in ten clads, including CYP51/72/74/85/86/97/710/711/727/746. 10 Gene duplications and divergence are responsible for this great diversity of CYP gene families under adaptive diversification. 11,12 In-plant genome, CYPs enzymes are more abundant compared to other organisms. Genome-scale identification of CYP has been performed in numerous model plants, such as 272 members in Arabidopsis thaliana, 312 members in Populus trichocarpa, 316 members in Vitis vinifera, 356 members in Oryza sativa, 372 members in Sorghum bicolor, and Glycine max 322 CYPs. [13][14][15][16] However, the functions of most CYPs are still uncertain. For instance, 30% of CYPs genes are well characterized in Arabidopsis. 17 CYPs play a pivotal role in lipid biosynthesis, fungal mycotoxins, and flavonoid biosynthesis. 18,19 A legume-specific CYP 93 C (CYP93C) is a key enzyme of isoflavonoid biosynthesis pathway. 20 Several cytochromes such as CYP93B/E/G play critical roles as plant secondary metabolites like triterpenoid saponin and flavone biosynthesis. 21,22 Similarly, CYPs have already been reported in plant organs such as germinating seed, stem, apical bud, and tulip bulbs. [23][24][25] 5-methylcytosine DNA methylation (5mC) is the most common epigenetic modifier and regulator of their targeted genes. In plants, 5mC and histone modification played vital roles in regulating a variety of physiological processes such as gene expression, genome stability, regulating parental imprinting, and in-plant response to abiotic or biotic stresses. 1,[26][27][28][29] However, the global DNA hypomethylation or hypermethylation was also estimated, but their exact mechanism remains unclear. [30][31][32] Rice (O. sativa L.) is a staple food crop that feeds half the world's population and is an ideal model plant owing to high sensitivity to salinity and drought, particularly at critical growth stages such as pollen development and seedling stage. 33,34 Rice plants can cope with environmental hazards by modulating various physiological and biochemical processes. These changes include the changes in enzymatic activities and expression level of different stress-responsive genes involved in salinity, drought, hypoxia, and temperature tolerance. 35 However, the functional characterization of these genes on behalf of DNA methylation yet needs to be explored.
Previously, we reported the dynamic role of DNA methylation 36 in regulating stressresponsive zinc-finger protein (ZFP) in rice and found that these genes have a varying degree of methylation. 37 Wei et al. 38 reported that the genome-scale identification of CYPs, but their role in rice DNA methylation regulation is still scarce. To avail this opportunity, we identified a total of 14 DNA methylation regulated CYPs in rice. It was observed that these genes are full or partially methylated in different regions and have the potential to tolerate abiotic stresses such as salt, drought, cold, and phytohormones in rice.

Phylogeny of CYPs in Rice and Arabidopsis
The phylogenetic relationships among the rice and Arabidopsis CYPs were analyzed in MEGA7 to generate an unrooted neighbor-joining (NJ) phylogenetic tree described by Wei et al., 38 Arabidopsis CYP protein sequences were retrieved from TAIR (http://www.arabidopsis.org/index.jsp) and aligned in ClustalOmega (https://www.ebi.ac.uk/Tools/ msa/clustalo/) and finally submitted in MEGA7. 41 The bootstrap was set at 1000 replicates.

Plant Growth Conditions and Material Collection
The surface-sterilized seeds of Japonica rice (O. sativa cv. Zhonghua11) were grown under conditions described in Ahmad et al. 37 The fresh samples of different plant parts, such as stem, root, leaf, node, and panicle, were harvested. For expression analysis of CYPs under various stresses, such as drought, salinity, and phytohormones, including indole-3-acetic acid (IAA), abscisic acid (ABA), and gibberellin (GA), two-week-old rice seedlings were used following as described in Ahmad et al. 37 Samples were collected at 0 h, 1 h, 3 h, 6 h, 12 h, and 24 h after post-treatments.

cDNA Synthesis and Quantitative Real-time PCR (Qrt-PCR) Assay
Trizol reagent (Invitrogen, USA) was used for total RNA extraction from harvested samples following the manufacturer's guidelines. DNase1 (Takara, Japan) was used to remove genomic DNA contamination from samples. RNA was further purified by precipitating the digested DNA. About 2 μg of total RNA was used for cDNA synthesis in the RT-PCR system (Promega, USA) following the manufacturer's instructions. To perform qRT-PCR of selected genes, SYBR Green Master Mix (Toyobo, Japan) was used in ABI-7500 Fast Real-Time PCR (Applied Biosystems, USA). QuantPrime tool 42 was used to design gene-specific primers (Table S1). However, 18s-rRNA was used to normalized gene expression. The expression was estimated by adopting ∆Ct method as described in Jain et al. 43

Statistical Analysis
GraphPad prism was used for statistical analysis. Three biological replicates for each sample were used for RT-qPCR analysis and five technical replicates were analyzed for each biological replicate. To determine significant differences two-tailed Student's t-tests was used at P values (*P < .05; ** P < .01; *** P < .001; **** P < .0001).

Identification and Characterization of Methylation Regulated Candidate CYP Genes
We adopted various bioinformatics approaches to search OsCYP of the grapevine genome, including BLAST and the hidden Markov model (HMM) profile method. A total of 384 CYP genes, including those 355 CYP genes previously identified by Wei et al., 38 were detected and subsequently verified for the P450 domain (PF00067). Among them, 14 potential CYPs showed varying methylation levels were identified through methylated DNA immunoprecipitation (MeDIP)-chip hybridization in rice under salt and osmotic stress. The protein length of methylation regulated OsCYPs ranging from 106 aa (LOC_Os09g26950) to 551 aa (LOC_Os08g43390). Table 1 shows the additional information about OsCYPs including their peptide length, pI, MW, insilico subcellular location, and genome distribution. According to gene mapping, 14 OsCYPs were anchored on six chromosomes in the rice genome ( Fig. 1). Chromosomes 1 and 6 anchored four genes each, three on chromosome 2, and a single gene was found on chromosomes 7, 8, and 9. Moreover, the number of exons varies from one to nine among 14 methylations regulated OsCYPs. For example, the single exon is found in LOC_Os09g26950 and LOC_Os02g12680. However, LOC_Os07g29960 and LOC_Os01g43740 displayed nine and five exons each whereas, the rest of the genes contained two exons except for LOC_Os02g17760, which contained three ( Figure S1).

Phylogenetic and Motif Analysis of Rice CYPs
To determine the phylogenetic relationship of OsCYP genes with CYPs from Arabidopsis, Solanum tuberosum, and P. trichocarpa, an unrooted NJ phylogenetic tree was generating. All CYPs were clustered in 9 clades (Fig. 2) and rice methylation regulated CYPs grouped in 6 different clans in the phylogenetic tree along with other OsCYPs (Fig. 2). Eight OsCYPs were clustered in Group-I, followed by two OsCYPS in cluster IV. Moreover, group II, VI, VII, and VIII contained a single rice protein. MEME motifs analysis revealed five unique motifs in 12 OsCYPs but absent in LOC_Os01g50520 and LOC_Os09g26950, while LOC_Os01g43720 contained three motifs only. LOC_Os02g12680 and LOC_Os01g43740 contained motif 2/4 and motif 1/3/4, respectively. LOC_Os07g29960, LOC_Os06g43440, LOC_Os06 g42610, and LOC_Os01g50530 showed motif 1/2/3/ 4 in their protein sequences. The rest of the OsCYPs contained all motifs in their peptides ( Figure S2).

DNA Methylation Status and Correlation Analysis of Candidate CYPs in Rice
By analyzing microarray data in three different conditions, 14 unique P450 genes exhibiting coordinately change in DNA methylation and expression were identified. To analyze whether the genes significantly affected epigenetically (>2-fold changes in methylation and expression) or not, a variable trend of methylation under two comparative conditions (control vs. stressed) was assayed. Here, an average normalized depth of reads as the measurement of methylation level for each gene was used. To quantify methylation level in control, salt, and osmotically stressed samples, we used MEDIPS software package. 44 To quantify DNA methylation in genes, all the CYPs were divided into different basic genic regions like a promoter, gene body, 5′ UTR, 3′ UTR, and downstream. Genome-scale peak scanning (unique peaks) for methylation was based on behalf of a defined analysis model (MACS), and the peak standard was fixed on behalf of P-value (p < 10e-5). Peaks of two samples (control vs. stressed) were merged as candidate different methylation regions (DMRs). Genomic regions with at least three consecutive windows statistically significant differentially methylated samples were considered differentially methylated regions (DMRs). Overall, we observed a wide range of methylation differences in CYPs under two comparative conditions, and the other was predicted in the promoters of the basic region of methylation, and out of 14 genes, nine were methylated in the promoter region. After promoter, downstream and gene body (CDS) were highly methylated regions among CYPs. Different quantities and varying levels of methylation in different regions represent the patterns of methylation, and the changed level of methylation might be the cause for their up/ downregulation. Among CYPs genes, we also observed the enrichment of methylation in different genes under two different stresses (control vs. salt or PEG6000 stress). LOC_Os02g12680 showed the methylation quantity in upstream 2 K and downstream 2 K region and showed downregulation under both stresses. Similarly, LOC_Os06g42610 showed the enrichment of methylation in downstream 2 K, CDS, and intron regions and downregulated under salt and osmotic stresses. Another gene, LOC_Os09g26950, was methylated in upstream 2 K and CDs regions.
Two different natures of stresses cause different changes epigenetically in the rice genome through the up/downregulation of genes. For example, the gene LOC_Os08g43390 showed its maximum DNA methylation 530 under salt stress compared to osmotic stress, where its methylation quantity was 444. Similarly, another gene, LOC_Os06g42610, showed 1125 quantities of methylation under salt stress and 1052 under osmotic stress, while its level of methylation was only 457 under control for downregulation. Another CYP gene (LOC_Os02g12680) was less methylated (98) under osmotic stress compared to salt stress (119). The detailed information of CYP genes about their methylation status and correlation under two comparative conditions is given in Table 2 and Fig. 3a-b.

Expression Analysis of OsCYPs in Different Parts/ organs of the Rice Plant
The specific expression pattern of methylationregulated OsCYPs was analyzed by qRT-PCR in different plant parts including, stem, root, node, panicle, and leaf. It was observed that OsCYPs are highly specific in their expression pattern (Fig. 4) which support the idea of tissue-specific epigenetic changes and the regulation of different CYPs under abiotic stresses. LOC_Os06g45960, LOC_Os02g12680, LOC_Os06g 43370, LOC_Os02g32770, LOC_Os01g43740, LOC_ Os06g42610, LOC_Os01g50530, and LOC_Os07g299 660 were highly expressed in root than other tissues. LOC_Os08g43390 expression peaked in nodes, while LOC_Os09g26950 have a weak expression in node but showed strong leaf expression. Similarly, LOC_Os01g50520 and LOC_Os06g43440 peaked in leaf, but the former has a weak expression in root. LOC_Os02g17760 and LOC_Os01g43720 were highly expressed in a panicle (Fig. 4).

Phytohormone Induced Expression Patterns of OsCYPs
Numerous cis-regulatory sequences related to hormone stresses were identified in promoter sequences of OsCYPs (Table S2). Thus, we examined the phytohormones-induced expression profile of these genes in the rice plant, including root and shoot (Fig. 5). The hormone-induced expression profile analysis revealed that few methylation-regulated rice CYPs genes were highly sensitive, while some of them were insensitive. LOC_Os06g43370, LOC_Os06g45960, LOC_Os07g29960, and LOC_Os08g43390 showed temporal expression at  all time intervals. It was found that the promoters of these genes had cis-regulatory elements related to ABA response (ABRE) (Table  S2). LOC_Os01g43740, LOC_Os02g32770, LOC _Os09g26950, and LOC_Os06g42610 showed peaked expression at three hours. However, LOC_Os08g43390 at six h; LOC_Os01g43720 and LOC_Os01g50520 at 12 h; LOC_Os02g17760 at 24 h were upregulated (Fig. 5a). For IAA treatment, LOC_Os01g50520, LOC_Os08g43390, and LOC_Os06g45960 were upregulated at three hours but downregulated during later time intervals. LOC_Os01g50530, LOC_Os06g43370, LOC_Os0 7g29960, LOC_Os09g26950, and LOC_Os06g42610 peaked at 12 h after IAA exposure. Similarly, LOC_Os01g43740 and LOC_Os02g32770 were upregulated at 24 h after treatment (Fig. 5b). Similarly, we observe that only of these genes' promoter contained cis-regulatory elements related to IAA among IAA-activated genes. These genes include LOC_Os02g32770 LOC_Os06g42610, LOC_Os06g45960, LOC_Os07g29960, and LOC_Os08g43390 (Table S2). For GA treatment, LOC_Os01g43720, LOC_Os07g29960, and LOC_ Os06g45960 were downregulated after treatment. LOC_Os01g43740, LOC_Os08g43390, LOC_Os01 g50520, and LOC_Os02g32770 were upregulated all time intervals except for LOC_Os08g43390 and LOC_Os01g50520, which was downregulated at 24 h after treatment. It is worth mentioning that only the promoter of LOC_Os01g50520 contained gibberellin-responsive elements (Table S2), suggesting that other regulatory genes may interact with these CYPs to activate their response to GA. However, LOC_Os06g42610 and LOC_Os09g26950 peaked at 24 h intervals (Fig. 5c). It is interesting to note that few genes were upregulated under all stresses while others induced under one of three stresses such as LOC_Os02g32770, LOC_Os06g42610, LOC_Os06g4337, and LOC_Os08g43390 was induced under all stresses whereas LOC_Os 01g43720 was induced under IAA, but LOC_Os01g43740 was suppressed. The phytohormones-induced expression of OsCYPs suggested that these stresses might regulate their expression pattern.

Expression Analysis of Rice Methylation Regulated CYPs Genes under Salt and PEG Stresses
To study the transcriptional regulation of OsCYPs under salinity and drought, we also analyzed the cis-acting elements in the promoters of CYPs. It was predicted that CYP genes contain cis-acting motifs such as temperature (LTR), drought (MBS), light (G-Box), and the pathogens may respond to abiotic stresses (Table S2). To understand the salinity and drought-induced expression profile of OsCYPs, the expression of 14 methylation-regulated genes was profiled. It was observed that LOC_Os01g43740 and LOC_Os09g26950 induced at 3 h; LOC_Os02g17760 and LOC_Os07g29960 at 6 h; LOC_Os06g42610 at 12 h; LOC_Os01g50530 at 24 h under salt stress.  h. LOC_Os02g32770 and LOC_Os09g26950 upregulated at 3 h; however, LOC_Os07g29960 was peaked at 6 h; LOC_Os01g50520 and LOC_Os06g43370 at 12 h; LOC_Os02g17760 was significantly peaked at 24 h (Fig. 6b). Overall, some genes were upregulated or downregulated under both stresses, while few showed opposite trends. LOC_Os01g43720 and LOC_Os06g45960 were downregulated under salinity, and LOC_Os02g17760, LOC_Os06g42610, LOC_Os06 g43370, and LOC_Os06g45960 were suppressed under drought. LOC_Os01g50530, LOC_Os06g4 3440, and LOC_Os06g43370 were downregulated under both stresses.

Discussion
DNA methylation (5mC) plays a dynamic role in regulating plant development, helping plant adaptation to environmental changes. [45][46][47] The mechanism by which DNA-methylation participates in these complex processes and mechanism is still unclear. The salinity and drought demethylate their responsive transcription factors, such as WRKY, NAC, TCP, and ARF, have been reported as key transcriptional regulators under stress conditions on behalf of DNA-methylation. [48][49][50] Nicotiana tabacum DNA-methyltransferase 1 (met1) mutant removes methyl-group from various genomics regions, activating stress-responsive genes. 51 However, the general mechanism of DNAmethylation associated with these stress-responsive genes remains unknown. Genome-scale DNAmethylation analysis was performed in various organisms from insects, chordates, animals to plants. 52,53 This revealed that it might be considered as the highlighted feature of all eukaryotes. 54 Rice is an important cereal crop worldwide, and a facile genomic model provides an opportunity to study the molecular biology of cereal plants. Rice DNA-methylation is of great interest to geneticists, biotechnologists, and breeders. We observed a significant enrichment of upregulated and downregulated genes associated with DNA-methylation or demethylation in rice. It was also analyzed that a large fraction of DEGs (differentially expressed genes) did not exhibit a significant difference in their methylation level (Table 2). Different researchers reported similar findings describing the role of DNA-methylation in transcriptional regulation of specific loci in different crops. 51,55,56 Genome analysis of these genes indicated that DNA-methylation plays a pivotal role in regulating hormone-specific CYP genes. A strong correlation of gene promotor and body methylation/demethylation for genes up/downregulation compared to any other flanking sequence is an essential feature of stress-responsive genes like CYPs to support the plant the adaptation of environmental stress like salt and drought. These findings are consistent with previous reports, including rice. 6,48,57,58 In the current study, the DNA-methylation in CYPs is estimated for the first time through genome analysis. Changes in the DNA-methylation level of CYPs and other related regulatory genes are the key to future crop breeding programs. Rice, being a saltsensitive cereal crop, displayed high genotypic variability. 59,60 These genotypes can tolerate extreme salinity, while few are highly sensitive. These are indicating that the unique rice genetics and regulatory mechanism are essential for further improvement. This study identified a unique array of genes, including stress tolerance, chromatin modification, retrotransposons, and hormone responsive. It suggests that the evolution of natural variability in rice germplasm for stress tolerance might be independent of the patterns or extent of DNA-methylation. Wang et al. 47 unveiled that DNA-methylation variation exhibited tissue specificity and caused significant differential gene expression. In tissues-specific expression analysis, we also find a diverse and variable expression of P450 genes in rice tissues that all genes are constitutively expressed (Fig. 4). Finally, to demonstrate transcription regulation, the cis-acting sequences analysis revealed that OsCYP genes respond to various abiotic stresses such as abscisic acid, MeJA, and ethylene lowtemperature, gibberellin, salicylic acid, drought, auxin (Table S2).

Conclusion
Conclusively, a dynamic correlation was observed between DNA methylation and CYP genes' expression under abiotic stresses by high throughput analysis. Varying level of DNAmethylation was observed in different regions of concerned genes under two comparative conditions. This may show their regulatory role in rice plants against these stresses and predicts the epigenetic role of DNAmethylation in gene suppression or activation. This study provides the foundation to explore the enzymatic activity of CYPs and find out the impact of DNA-methylation in gene stress adaptation mechanisms for rice and other related crops to improve growth and yield as well.