Genome-wide identification and expression profiling reveal the diverse role of Methyl-CpG-binding domain proteins in tea plant Camellia sinensis

As a conserved epigenetic mark, DNA methylation is dynamically regulated by DNA methylases and demethylases, thereby playing critical roles in gene regulation^[1], transposon silencing^[2], and chromosome interactions^[3]. In plant genomes, DNA methylation is more likely to occur in three types of bases (CG, CHG and CHH) during development and reproduction^[4]. In plants, symmetric CG and CHG methylation occurs during DNA replication, which is catalysed by DNA METHYLTRANSFERASE 1 (MET1) and CHROMOMETHYLASE 3 (CMT3), respectively^[4]. In contrast, asymmetric CHH methylation is established through de novo methylation by DOMAINS REARRANGED METHYLTRANSFERASE 2 (DRM2)^[5]. However, DNA glycosylase/lyases can catalyse DNA demethylation, with REPRESSOR OF SILENCING 1 (ROS1) and DEMETER (DME) being typical representatives^[4]. Methyl-CpG-binding domain (MBD) proteins recognise DNA methylation signals and bind methylated DNA to contribute to transcriptional repression, chromatin remodelling, and DNA demethylation^[6−8]. The first plant MBD protein was isolated from pea (Pisum sativum)^[9]. Subsequently, the MBD family members were identified in Arabidopsis, rice, maize, tomato, and wheat^[10−12].

In Arabidopsis thaliana, 13 members from the MBD family were classified into eight groups according to the similar MBD motif^[7]. The function of most members is unclear, but three (MBD5, MBD6, and MBD7) of them were demonstrated to have capacity for binding symmetrically methylated DNA^[7]. For example, AtMBD6 is the interactor of RNA binding proteins and histone deacetylase (AtHDA6, histone deacetylase 6) which regulates the RNA-mediated gene silencing^[13]. AtMBD7 can form the MBD7-IDM complex with IDM1/2/3 (histone acetyltransferase IDM1, alpha-crystallin domain proteins IDM2 and IDM3) and HDP1/2 to regulate the recruitment of ROS1 promoting DNA demethylation^[14]. In maize, 14 MBD genes have been identified^[12], including ZmMBD101, which were reported to have capacity for binding DNA and maintaining mutator elements of chromatin in an inhabited state^{[15, 16]}. Moreover, in maize, the ZmMBD11 also showed specific DNA-binding activity indicating its important role in reading cytosine methylation^[12].

Rice contains 17 proteins with MBD domains, in which the function of partial genes has been reported. The deposit of cyclobutane pyrimidine dimers were reduced by OsMeCP (MBD protein) on the rice genome^[17]. While in tomato, MBD proteins were found to participate in fruit ripening^[11]. These studies in different plants indicate a significant role for MBD proteins in genomic and genetic regulation. However, to date, the potential role of each individual member of the MBD family has not been demonstrated in any perennial woody plants.

The tea plant (Camellia sinensis (L.) O. Kuntze) is a perennial woody cash crop with a wide distribution area. During its growth cycle, C. sinensis encounters biotic and abiotic stresses including cold, drought, plant diseases and pests. During cold acclimation, a change in the DNA methylation level on the tea plant genome was suggested to be involved in cold resistance in this species^[18]. The relative transcription levels of genes encoding DNA methylase and demethylase showed significant change in response to abiotic stress^[19]. Epigenetic regulation was identified as one of the major mechanisms involved in overwintering bud dormancy formation and release^[20]. However, the specific factors in the epigenetic regulation system in tea plants is largely unknown, especially the response of DNA methylation to adverse environmental conditions, growth, and development. Therefore, we present an overall study of DNA methylation reader proteins in tea plants. Since the draft genome of 'Yunkang 10' was reported, a total of eight genomes of different tea plant species have been sequenced and released. Among them, the 'Longjing 43' (LJ43) genome was one of the high-quality chromosome-scale tea genome, which was characterized by a scaffold N50 value of 144 Mb, 88.36% gene completeness, and a base accuracy of 99.999%. It provides a rich genomic resource for tea researchers^[21]. Thus, we identified and characterised 16 genes encoding MBD proteins within the LJ43 genome^[21]. Comprehensive expression detections of these genes under environment stresses and different growth and development stages were performed to provide useful insights into the epigenetic regulation in tea plants.

The Hidden Markov Model (HMM) profile of the MBD domain (PF01429; http://pfam.sanger.ac.uk) was used to scan predicted proteins encoded in the MBD genes in the LJ43 genome (https://bigd.big.ac.cn/search/?dbId=gwh&q=GWHACFB00000000^[21]). All candidate proteins were determined using the SMART database (https://smart.embl-heidelberg.de). The potential members of the CsMBD gene family were confirmed by homology analysis using the BLAST algorithm (p-value < 10⁻⁵) combined with MBD domain identification. The newly identified genes were named according to the Arabidopsis gene locations on the phylogenetic tree. The WoLF PSORT server (https://wolfpsort.hgc.jp) was used to predict the subcellular localisations of these proteins. The physico-chemical parameters of these gene encoding proteins were calculated using the ProtParam tool (https://www.expasy.org).

Phylogenetic tree construction and analysis of sequence alignment, conserved domains, gene structures, chromosomal location, duplication analysis and regulatory elements

To further determine the phylogenetic relationship of MBD proteins, the conserved MBD sequences from four plants (Oryza sativa, A. thaliana, Populus trichocarpa, C. sinensis) were aligned using ClustalW^[22]. The unrooted tree was created using MEGA7.0 (https://www.megasoftware.net) by the neighbour-joining (N J) method. The bootstrap value = 1,000. The conserved domains of these proteins were analysed using Tbtools software^[23] and the chromosomal location of the genes was investigated and analysed by MG2C (http://mg2c.iask.in/mg2c%5Fv2.1) based on their starting positions on the chromosomes. Exon/intron structure of CsMBDs were predicted with the Gene Structure Display Server 2.0 (http://gsds.gao-lab.org)^[24]. The genome sequences and gene annotation files were used to conduct the duplication by TB tools. The 2 kb long upstream regions were chosen to analyse cis-acting elements using the Plant CARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html).

Plant materials and treatments

LJ43 cultured at the Tea Research Institute of the Chinese Academy of Agricultural Sciences (Hangzhou, Zhejiang Province, China) was selected. To examine tissue-specific gene expression patterns, the leaves, roots, buds, and stems of LJ43 were collected in August 2020. Dormant axillary buds in the middle of branches were sampled from October 2018 to March 2019. In this study, the paradormant axillary buds were collected on 25 October 2017. The endodormant axillary buds were collected on 14 January 2018 and the ecodormant axillary buds were collected on 13 March 2018. Cold acclimation leaf sampling included mature LJ43 leaves sampled on 1 November and 15 November 2017, 18 December 2017, 11 January 2018, 6 February 2018, and 27 March 2018. All experimental materials and treatment are shown in Supplemental Table S1.

Five-year-old LJ43 cutting seedlings were planted in a climate chamber at 25 °C in a 16/8 h light/dark cycle with 60% humidity. For drought treatment, water was provided every 3 days to the control group, while the experimental group was given no water. The second and third leaves of the control and drought treatment plants were sampled on 0, 3, 6, 9, 12, 15 and 18 d. For seeding germination assays, the seeds were stripped of the seed coat and grown on vermiculite before being uprooted after 2 weeks. For floral development analysis, flower buds, semi-open flowers, and fully open flowers were collected. For biotic stress treatment, healthy branches were inoculated with conidial suspensions of Colletotrichum camelliae (10⁶ spores/mL) as previously described^[25]. The inoculation method for Ectropis oblique and Pseudopestalotiopsis camelliae-sinensis was conducted as previously described^[26], while the control groups were given a sterile water spray. The bud and first leaves were sampled at 12 and 24 h post-inoculation.

Three biological replicates were set for the above sampling and treatment. After sampling, the samples were immediately placed in liquid nitrogen and stored at −80 °C for later use.

Gene expression level detection

Total RNA was isolated from the experimental samples using an RNAprep PurePlant Kit (Tiangen, Beijing, China) by following the manufacturer's instructions. All cDNA samples were synthesised using PrimeScript RT enzymes with the gDNA eraser (Takara, Dalian, China). The qRT-PCR reactions were set up using LightCycler® 480 SYBR Green I Master (Roche Diagnostics GmbH, Mannheim, Germany). All primer pairs for expression detection are shown in Supplemental Table S2 and CsPTB1 was chosen as the internal reference gene^[27]. The $2^{-\Delta} $^Ct method was used to calculate the expression levels of cold acclimation‒related and tissue-specific genes, while the $2^{-\Delta \Delta} $^Ct method was employed in the relative expression analysis of the remaining experimental treatments^[28]. Two technical replicates were performed for each biological replicate.

Cloning for full-length open reading frames and promoter sequence validation

Total RNA of LJ43 buds was extracted as described previously. The synthesise of cDNA is performed using the DNase I Amplification Grade and SuperScript^@ III First-strand kit (Invitrogen, USA) following the manufacturer’s instructions. Based on previous data analysis for the LJ43 tea genome^[21], the full-length open reading frames (ORF) encoding the MBD genes in tea plants were amplified using the primers provided in Supplemental Table S3. PCR amplification was conducted in 50 μL reaction mixtures for KOD-plus-Neo (TOYOBO, Osaka, Japan). The PCR product was gel-purified and cloned into the pEASY-Blunt zero vector (TransGen Biotech, Beijing, China) for sequencing by the Ykang Biotechnology Company (Ykang, Hangzhou, Zhejiang). The obtained sequences were uploaded to the National Center for Biotechnology Information (NCBI).

The genomic DNA of the LJ43 buds was isolated using the Plant Genomic DNA Kit (Tiangen, China). The upstream regulatory elements present in the first 2,000 bp of the start codon were used as the gene promoter for cloning. The primers are listed in Supplemental Table S4. KOD-plus-Neo (TOYOBO, Osaka, Japan) was used in PCR amplification. The PCR products were gel-purified and cloned into the pEASY-Blunt zero vector (TransGen Biotech, Beijing, China) for sequencing by the Ykang Biotechnology Company (Ykang, Hangzhou, China). The cis-elements in the obtained sequences were identified and uploaded to the NCBI (accession numbers shown in Supplemental Table S4).

Data statistical analysis and graphic production

TTEST function in Microsoft Excel XP was used for statistical differences analysis, where p < 0.05 was considered significant and p < 0.01 was considered highly significant. Values are expressed as mean ± standard deviation. Except for the heatmaps by TBtools^[23], all graphs were generated using GraphPad Prism 8 (https://www.graphpad.com).

In this study, 16 genes encoding MBD domain-containing proteins were identified from the LJ43 genome databases. The basic characteristic analysis showed that the ORF of 16 CsMBD genes varied in size. The protein length of CsMBDs varied from 128 to 1,317 residues, with an average length of approximately 359 amino acids. Among these proteins, the smallest MBD was 128 amino acids long (CsMBD6), while the largest protein was composed of 1,317 (CsMBD9). The predicted molecular weight for CsMBD proteins ranged between 14.31 and 142.71 kDa, with an average molecular weight of approximately 47.58 kDa. The theoretical PI was between 4.54 and 9.29, while the theoretical PI values of nine proteins were smaller than 7.0. With the exception of CsMBD8, the instability index for MBD proteins isolated from tea plants was more than 40, indicating that these predicted proteins were unstable. Additionally, grand average hydropathicity (GRAVY) analysis showed that all CsMBD proteins were hydrophilic and according to their subcellular localisation, they were all localised in the nucleus (Table 1).

For cloning validation, we cloned the ORF of MBD gene family, with the exception of CsMBD7, we only obtained its partial sequences from genome DNA. Moreover, all the obtained sequences have been uploaded to the NCBI (accession number as shown in Table 1).

Analysis of tea plant MBD protein structures

The conserved domains of MBD proteins in tea plants and Arabidopsis were determined using the SMART database and full protein sequences. According to structural domain analysis, 14 different motifs were determined in CsMBD proteins, while all CsMBD proteins contained the MBD domain (PF01429). Sixteen MBD proteins in tea plants were divided into eight classes according to the differences in the MBD domain in the sequences (Fig. 1) . All MBD proteins from class II and III contained a zf-CW domain related to transcriptional cis-acting element binding. CsMBD9 was placed in class V, which includes proteins with the most motifs, including the PHD, MBD, DDT and HSP70 domains. Interestingly, AtMBD13 fell into class VIII, which includes proteins that lack some domains. However, CsMBD13 and CsMBD14 were placed in this class even though they contain several domains. Moreover, CsMBD15 was most similar to AtMBD7, because they have two MBD motifs that share a high sequence similarity, suggesting functional equivalence.

Figure 1.  Schematic structures of MBD proteins; CsMBD1-16 from the tea plant and AtMBD1-13 from Arabidopsis. The conserved MBD domains are indicated with green rectangles. The length of the box represents the length of the domain. Different coloured rectangles represent different protein domains.

To study the structural characteristics of the MBD proteins, the conserved MBD sequences from three species (tea plant, human, Arabidopsis) were aligned (Fig. 2). The results showed a high homology of amino acid sequences between the arabidopsis MBD motif of AtMBD1–AtMBD13 and that of the tea plant CsMBD1–CsMBD16. Here, we focused on these residues in the methyl cytosine binding and guanine recognition sites. Five residues (Arg22, Arg30, Asp32, Tyr34 and Arg44) were reported to be important for MBD protein DNA binding^[29]. As shown in Fig. 2, these five residues were largely confined to the MBD motif both in the Arabidopsis and C. sinensis proteins. Indeed, the MBD protein sequences in the same class were also highly similar. Moreover, MBD proteins consist of four β-sheets, an α-helix, and two hairpin loops, which form the MBD–C^mDNA complex. The recognition for the methyl group depends on conserved residues within the β-sheets, while the other structures are in contact with the DNA backbone^[29].

Figure 2.  Alignment of MBD domains of the AtMBD1-13 from Arabidopsis, CsMBD1-16 from the tea plant and MeCP2 proteins from humans. The secondary structure for amino acid conservation is displayed at the top: four β-sheet strands (β1–4), one α-helix, and two loop regions (L1–2). Five asterisks indicate the positions of residues that may be crucial for the binding of MBD proteins to methylated DNA.

Phylogenetic analysis and duplication analysis of CsMBD genes

To further understand the evolutionary relationships of MBD proteins, related MBD domain sequences from four species (O. sativa, Populus, C. sinensis and A. thaliana) were chosen to generate the phylogenetic tree (Fig. 3). The result showed that the MBD proteins of these four species can be classified into eight classes. Class I and VI contained the most MBD proteins from C. sinensis, while class IV and VIII were comprised of two C. sinensis proteins. The other classes contained one protein from the tea MBD family. Phylogenetic analysis revealed sister pairs of paralog genes; for example, one At-Cs sister pair (AtMBD10 and CsMBD10), one At-Pt sister pair, and one Pt-Cs sister pair. However, sister pairs from orthologous genes were more common than those from paralogous genes, which indicated that the orthologous genes may have descended from a common ancestor.

Figure 3.  Phylogenetic tree of MBD proteins in tea plants (16), poplar (15), rice (14) and Arabidopsis (13). The unrooted tree was created using MEGA7.0 by the neighbour-joining method. The bootstrap value = 1,000.

To investigate the genomic distribution of the MBD gene family in C. sinensis, the physical chromosomal locations of CsMBD genes were revealed. Fifteen MBD genes were distributed unevenly across nine of the 15 chromosomes (Fig. 4a). Notably, one gene (CsMBD2) was not mapped within the chromosome while it was assembled to the scaffold. Chromosome 11 contained the maximum number of genes (four), which were located in the middle of the chromosome. Chromosomes one, three, four, six, and eight only contained one MBD gene each, which were distributed randomly on the chromosomes. Chromosomes 12, 13, and 14 each had two MBD genes that were located in the middle and end of the chromosomes. Moreover, nine paralogous pairs were found in the CsMBD family genes (Fig. 4b). For example, CsMBD11 and CsMBD12 form a pair of orthologous pair. The syntenic analysis of MBD genes showed that 16 CsMBD genes were collinear with 8 Arabidopsis genes.

Figure 4.  Duplication analysis of MBD genes in tea plants and Arabidopsis. (a) Synteny analysis of 16 CsMBD genes in tea plants. (b) Synteny analysis of MBD genes between tea plants and Arabidopsis. Gray lines in the background indicate the collinear blocks within the tea plant and other plant genomes, while the red lines highlight the syntenic MBD gene pairs.

Cis-acting element analysis of CsMBD genes

To study the potential biological functions of CsMBD genes, regulatory elements were identified in the promoter region with a length of 2,000 bp (Fig. 5). The cis-acting elements of CsMBD genes can be classified as light, phytohormone, stress responsiveness, and plant growth and development-related elements (Fig. 5). Among these four categories, the light responsiveness classification contained the most cis-acting elements, including an ACA-motif, ACE, and AE-box. With the exception of CsMBD6, CsMBD8, and CsMBD16, the other 13 genes contained four responsive classifications, among which CsMBD11 had the largest number of stress-responsive elements. However, CsMBD8 contained the most cis-acting elements without any regulatory elements for plant growth and development. In CsMBD6, we did not find any plant growth or development-related elements. In total, light and phytohormone responsive elements were the most common in all the CsMBD genes evaluated.

Figure 5.  Analysis of cis-acting elements in the promoter region of 16 CsMBD genes. (a) The number of different cis-acting regulatory elements in the promoter region. (b) The histogram of cis-acting elements in four categories.

Tissue-specific expression of CsMBD genes

To gain insight into their putative functions, we first analysed their expression patterns in different tissues (Fig. 6). The expression profiles of CsMBD genes in tissues or organs were mainly divided into three types according to cluster analysis. The first expression pattern included nine CsMBD genes, which had an extremely high expression in the buds and leaves. CsMBD2/8/1 had the highest expression in the buds, and CsMBD14/6 had the highest transcription in the leaves. The second type of expression pattern only included two genes (CsMBD13/9), whose expression was highest in the roots, followed by buds and leaves, and thus the lowest in the stems. CsMBD4/10 belongs to the third expression pattern, and the expression level was the highest in the stems, followed by roots, buds, and leaves. Additionally, the expression of CsMBD7 was not detected in the experimental samples. These results indicate that the CsMBD gene has three different expression patterns and obvious tissue specificity, and may play different biological roles in different tissues.

Figure 6.  Tissue-specific expression profile clustering of 15 CsMBDs. The colour scale indicates log₂ signal values (right); blue and red colours indicate low and high levels of transcript abundance, respectively. The different organs are noted on the right. A cluster dendrogram is shown above.

Expression changes of CsMBD genes in response to biotic stresses

After C. camelliae infection, ten genes out of 15 CsMBD genes (CsMBD1, CsMBD2, CsMBD4, CsMBD5, CsMBD6, CsMBD10, CsMBD11, CsMBD12, CsMBD15 and CsMBD16) were significantly upregulated after 12 h (Fig. 7) . The relative expression of other CsMBD genes showed no significant difference after 12 h. After Ectropis obliqua inoculation, the expression of nine genes (CsMBD1, CsMBD2, CsMBD3, CsMBD5, CsMBD9, CsMBD10, CsMBD13, CsMBD14 and CsMBD16) was significantly high after 12 h, while only three genes (CsMBD3, CsMBD5 and CsMBD9) were significantly upregulated after 24 h. However, except for CsMBD3, CsMBD4, CsMBD5, CsMBD8 and CsMBD10, we observed no significant changes in the transcript abundance of other genes between the infected and control plants after the Pseudopestalotiopsis camelliae-sinensis infection. Interestingly, the expression of CsMBD3 was significantly induced by Ectropis obliqua and Pseudopestalotiopsis camelliae-sinensis, while the expression of CsMBD8 was significantly induced by Pseudopestalotiopsis camelliae-sinensis infection.

Figure 7.  Expression profiles of 15 CsMBD genes under biotic stress treatments in tea plants. The different colours indicate different plant diseases and insect pests. The error bars indicate the standard error of the mean (SEM) from three biological replicates. * Significant differences (p < 0.05), ** highly significant differences (p < 0.01).

Expression profiles of CsMBD genes under abiotic stress conditions

Specifically, ten genes showed no significant changes in expression in response to cold stress, while five genes (CsMBD1, CsMBD5, CsMBD11, CsMBD15 and CsMBD16) were significantly induced (Fig. 8a). Under drought conditions, we observed that the genes of this family showed four expression patterns. The transcriptional abundance of six CsMBD genes (CsMBD1, CsMBD5, CsMBD11, CsMBD15, CsMBD13 and CsMBD9) reached its maximum after 18 d of drought (Fig. 8b). Moreover, the expression of five genes (CsMBD14, CsMBD2, and CsMBD10) were significantly upregulated after 15 d of drought. The transcript abundance of CsMBD6 and CsMBD12 were significantly upregulated after 12 d, while the transcript abundance of CsMBD4 and CsMBD8 showed no significant changes after being subject to drought conditions. Interestingly, CsMBD9 expression was upregulated by drought treatment after 9 d, with a relative expression level that was approximately 80 times higher than that measured at 0 d, and it remained high until 18 d. These findings revealed that most of CsMBD genes were induced by drought stress.

Figure 8.  Relative expression of 15 CsMBD genes under abiotic stresses in tea plants. (a) Cold acclimation; (b) drought stress conditions. The y-axis represents the transcript abundance, and the x-axis represents the stress treatments (days or date). The error bars indicate the standard error of the mean (SEM) from three independent biological replicates. * Significant differences (p < 0.05), ** highly significant differences (p < 0.01).

Expression analysis of CsMBD genes related to growth and development

Among the 15 genes evaluated, the expression profiles of three genes (CsMBD2, CsMBD3, and CsMBD4) showed a significant increase in expression during seed germination, whereas the relative expression of CsMBD10 decreased (Fig. 9). During floral development, the majority of CsMBD genes were downregulated, with only CsMBD9 showing no obvious difference in expression during the different stages of flower development (Fig. 9). Furthermore, CsMBD4 and CsMBD12 expression increased, while the other 12 genes decreased in expression. These findings revealed that the CsMBD genes exhibit different expression patterns during floral development.

Figure 9.  Expression patterns of 15 CsMBD genes during growth and development in tea plants. Including seed germination, floral development and different dormancy stages in tea plants; p < 0.05 was considered significant (*) and p < 0.01 was considered highly significant (**). The colour scale indicates the log₂ signal values shown on the right. Red and blue indicate high and low levels of expression, respectively. A cluster dendrogram is shown on the left of the graph.

During winter dormancy, the expression of CsMBD genes in axillary buds exhibited two different patterns (Fig. 9). The results showed that 6 CsMBD genes (CsMBD1/8/12/5/11/6) were expressed at their lowest levels in axillary buds during the endodormancy stages, while their expression levels were significantly increased during the ecodormancy stages. This expression pattern is highly correlated with the physiological state of bud dormancy. These results indicated that CsMBD genes are involved in bud dormancy in tea plants, especially in the process of bud dormancy release.

DNA methylation is essential in the response to environmental stress as well as to the developmental process. To better understand epigenetic regulation associated with environmental stimuli and plant growth and development, it is essential to study the methylation of MBD proteins. For now, only a limited number of MBD genes have been identified in plants, including 13 in Arabidopsis, 18 in tomato^[11], 14 in maize^[12], and 17 in rice^[11]. In our study, a total of 16 MBD genes were identified and characterised in the tea plant genome, which possess typical characteristics of MBD proteins. These identified proteins were divided into eight subclasses based on phylogenetic analysis, and were consistent with Arabidopsis and populus. Moreover, genes belonging to the same subclass share similar structures. For example, class IV, V, and VII, Arabidopsis and tea plants owned equal numbers of MBD genes. However, subclass V included an MBD protein found in dicots but not monocots^[10]. It is possible that either the MBD gene was inserted into the dicot genome, or lost from the monocot lineage. These results revealed that the evolution of MBD genes is related to the difference between monocots and dicots. In spite of the fact that plant genomes vary widely in size, such as the tea genome (approximately 3 Gb), rice (approximately 389 Mbp), and Arabidopsis (approximately 125 Mbp), the number of MBD genes were similar among these species, suggesting that MBD proteins are highly conserved throughout evolution. According to the synteny analysis of MBD genes between tea plants and Arabidopsis, we speculated some multicope genes and unicope genes during the evolution of tea plants. AtMBD1 and AtMBD3 were lost in the long process of evolution, which might be due to other DNA methylation reader genes taking over the biological functions of these two reader genes, or another alternative mechanism for DNA methylation. While four genes (CsMBD1, CsMBD4, CsMBD8 and CsMBD10) were novel genes, which suggests that these genes may have evolved new functions or a more complex mechanism.

Cis-acting element analysis defined four categories of elements in CsMBDs, including stress, phytohormone, light, and developmental response. These regulatory elements suggest that MBD genes have diverse roles in biological processes. Among these cis-acting elements, the TCA-element, involved in salicylic acid responsiveness, is commonly found in most MBD gene promoter regions and plays an indispensable role in response to biotic stress. In addition, RY-elements, related to seed-specific regulation, were found in CsMBD2 and CsMBD13 promoter regions, which indicates that these two genes are likely to be important in seed development. In tomato, the accumulation of SIMBD transcripts gradually decreased in wild types during fruit ripening and ethylene-responsive cis-elements were commonly found in these gene promoters, suggesting that these genes can be induced by exogenous ethylene^[30]. Furthermore, the overexpression of SlMBD5 caused pleotropic developmental changes, including darker green fruits due to the accumulation of plastids and elevated pigmentation, as well as an increased abundance of SlGLK2, which is essential in the regulation of plastid biogenesis^[31]. In Arabidopsis, the overexpression of AtMBD8 resulted in delayed flowering during both long and short days in the C24 background^[32]. Moreover, the atmbd9 mutant presented multiple phenotypes such as early flowering and the promotion of branches. These findings indicate that MBD genes participate in the regulation of plant growth, as well as in the responses to environmental stimuli.

Genes encoding MBD proteins are involved in biotic stress responses

Plants showed genome-wide DNA methylation alterations when they were infected with pathogens and colonised by microorganisms. In A. thaliana, mild but widespread differential DNA methylation was observed after exposure to pathogens^[33]. The DNA demethylase triple mutant ros1–dml2–dml3 was more susceptible to Fusarium oxysporum than the wild type A. thaliana^[34]. Moreover, most genes of biotic stress responses are repressed in ros1–dml2–dml3 plants, potentially because their promoters are enriched with small transposons^[34], which suggests that DNA methylation plays an essential role in regulating the transcription of defence genes^[5]. Among the biotic treatments tested, most CsMBD genes were significantly induced by Colletotrichum camelliae and were upregulated in the leaves 12 h after inoculation, indicating that these genes respond to C. camelliae infection by regulating DNA methylation. CsMBD3 and CsMBD8 are most significantly induced by a P. camelliae-sinensis infection; the promoter of CsMBD8 contains eight cis-acting elements involved in methyl jasmonate-responsiveness (MeJA), such as the GGTCA and TGACG motifs. MeJA can enhance the defensive compounds synthesis and initiate transcription of pathogenesis-related genes for resistance^[35]. Therefore, CsMBD8 is a candidate gene for understanding the epigenetic regulation system that is involved in biotic stress in tea plants.

Putative genes in the MBD family respond to abiotic stress

Numerous studies have shown that plants utilize epigenetic mechanisms to fine-tune their responses to abiotic stress^[36]. Chromatin interaction, histone modification, and DNA methylation during cold and drought stress have been studied extensively. For example, in A. thaliana, the chromatin remodelling protein, PICKLE (PKL), is involved in the CBF-dependent cold stress response^[37]. Moreover, the expression of two cold-responsive genes, COR15A and GALACTINOL SYNTHASE 3 (GOLS3), are accompanied by a decrease in H3K27me3 deposition in cold stress conditions^[38]. Histone acetylation also plays an important role in the cold response; HvDME was found to be induced in response to drought treatment in barley^[39]. We found that four genes (CsMBD1, CsMBD5, CsMBD11, and CsMBD16) were significantly induced by cold acclimation and most of the CsMBD genes were dynamic regulated between 0 and 18 d of drought stress, except for CsMBD1, CsMBD5, and CsMBD11. Interestingly, the transcript abundance of CsMBD5 was upregulated in response to both cold and drought stress. CsMBD5 contains three LTR and two ABRE cis-acting elements that are involved in low temperature responsiveness and abscisic acid responsiveness, respectively. Therefore, CsMBD5 most likely participates in cold and drought responses. Moreover, the expression level of all DNA demethylase genes in tea plants significantly increased at certain times between 0 and 48 h of drought stress and the expression of genes encoding DNA methyltransferase was significantly inhibited^[19]. These results demonstrate that DNA methylation may regulate the abiotic stress responses in tea plants. There is no doubt that DNA demethylase and genes encoding MBD proteins are differentially regulated. The inhibition of DNA methyltransferase promotes the demethylation of gene abiotic stress-related responses, induces their expression and promotes abiotic stress tolerance.

CsMBD genes are involved in tea plant growth and development

DNA methylation is essential for regulation in plant growth, but how MBD proteins participate in these processes is complex and the mechanism is still unclear. In Arabidopsis, transposable elements cause significant accumulation and loss of CHH methylation during seed development and germination, respectively^[40]. Those genes involved in RdDM and CMT2 pathways are highly expressed during seed development, and DNA demethylases are not expressed during germination^[40]. In tea plants, three CsMBD genes (CsMBD2, CsMBD3, and CsMBD4) were upregulated during seed germination. Typically, CsMBD2 is expressed only in buds and our promoter analysis revealed that the CsMBD2 promotor contains cis-regulatory elements related to endosperm expression and seed-specific regulation, such as the RY-element and the GCN4 motif. These results suggest that CsMBD2 may regulate seed germination through participation in the histone acetylation pathway. During Arabidopsis flower development, the CG, CHG, and CHH sites were methylated and accounted for 17.5%, 13.7% and 5.2% of the 2,4035 genes, respectively^[41]. Numerous cytosine sites were methylated during the meristem to the early flowering developmental stages^[41]. Consistent with our study, the expression of most MBD protein-encoding genes in tea plants was suppressed in the flower bud to the semi-open flower bud stages.

Epigenetic processes play crucial roles in the regulation of the transition between bud set and burst^{[42, 43]}. Genomic DNA methylation patterns are associated with the state of shoot apical cells when poplar shifts to the dormant stage from the growing stage, with PtaDML10 expression levels the highest during dormancy release^[43]. In tea plants, five CsMBD genes were maximally expressed in the bud break stage. Among these genes, CsMBD8 was specifically expressed only in the buds and was shown to contain eight ABRE cis-acting elements related to abscisic acid responsiveness (ABA). ABA is a widely reported exogenous hormone that regulates bud dormancy^[44]. The results indicate that MBD proteins may be engaged in the regulation of the growth cycle from winter dormancy to bud growth in plants.

The characterisation, features, and amino acid conservation of the CsMBD genes have provided new possibilities for improving light, phytohormone, stress, and plant growth and development responsiveness in tea plant crops. Comprehensive gene expression profiling of all CsMBD genes in tea plants was performed in different tissues and gene expression patterns were measured in leaves from plants grown under abiotic and biotic stress during different stages of development. We identified several CsMBD genes as candidates for further investigation due to their performance during plant stress responses and plant growth development. CsMBD8 could be a candidate gene for investigating the epigenetic mechanism of abiotic stress and bud dormancy, while CsMBD5 could serve as a model for understanding epigenetic mechanisms that are related to abiotic stress. In addition, CsMBD2 could also be used to study the epigenetic mechanisms related to seed germination. Overall, these findings provide valuable clues for identifying the specific functions of this gene family and the genomic diversity among different tea plant genotypes (Supplemental Fig. S1).

Data availability statement

The accession number of CsMBD genes sequences were obtained as shown in Supplemental Table S2, and accession number for cloning promoter sequences of CsMBD genes were exhibited in Supplemental Table S4.