Genome-wide identification and expression analysis of SABATH methyltransferases in tea plant (Camellia sinensis): insights into their roles in plant defense responses

ABSTRACT

SABATH methyltransferases convent plant small-molecule metabolites into volatile methyl esters, which play important roles in many biological processes and defense reactions in plants. In this study, a total of 32 SABATH genes were identified in the Camellia sinensis var. sinensis (CSS) genome, which were renamed CsSABATH1 to CsSABATH32. Genome location annotation suggested that tandem duplication was responsible for the expansion of SABATH genes in tea plant. Multiple sequence alignment and phylogenetic analysis showed that the CsSABATHs could be classified into three groups (I, II and III), which were also supported by gene structures and conserved motifs analysis. Group II contained only two CsSABATH proteins, which were closely related to PtIAMT, AtIAMT and OsIAMT. The group III SABATH genes of tea plant exhibited expansion on the CSS genome compared with Camellia sinensis var. assamica (CSA) genome. Based on RNA-seq data, the CsSABATHs exhibited tissue-specific expression patterns, and the members with high expression in buds and young leaves were also obviously upregulated after MeJA treatment. The expression of many transcription factors was significantly correlated with that of different members of the CsSABATH gene family, suggesting a potential regulatory relationship between them. Quantitative real-time PCR (qPCR) expression analysis showed that CsSABATHs could respond to exogenous JA, SA and MeSA treatments in tea plants. RNA-seq data analysis and qPCR validation suggested that CsSABATH8, 11, 16, 25, 29 and 32 might play a special role in plant defense against insect herbivory. These results provide references for evolutionary studies of the plant SABATH family and the exploration of the potential roles of CsSABATHs in tea plant defense responses.

KEYWORDS:

• Camellia sinensis
• SABATH gene family
• methyltransferase
• plant defense
• MeSA

Introduction

Tea plant [Camellia sinensis (L.) O. Kuntze] is an important perennial evergreen woody economic crop. This species is also sensitive to many fungi, bacteria, viral diseases and herbivorous pests, and such biotic stresses have a serious impact on the yield and quality of tea.¹ How to effectively prevent the loss of tea production caused by biological stress is an urgent problem to be solved. Over the long process of evolution, plants have formed complex environmental adaptation mechanisms, which can show themselves and their state to surrounding organisms through a variety of ways, the most important of which is the release of volatile substances.^2,3 Although plants can release volatiles before insect infestation, which plays a role in direct defense, the more widely occurring phenomenon in nature is that attack by pests can increase the release of volatiles and change their composition, which improves plant defensiveness as well as attracts the natural enemies of herbivores to play an indirect defense role.⁴ These kinds of volatiles are called herbivore-induced plant volatiles (HIPVs). The HIPVs can not only protect attacked plants from further damage but also act as warning signals to nearby plants, which in turn stimulate their defense responses.⁴ Therefore, HIPVs are effective for the biological control of tea plantation pests.

Volatile methyl esters, such as methyl salicylate (MeSA) and methyl benzoate (MeBA) are HIPVs.⁵ In plants, these compounds not only are the components of floral fragrance but also play a special role in plant defense.^6-8 They have been detected in many herbivore-attacked plants, including Arabidopsis thaliana, lima bean, tomato, soybean, rice, tobacco and tea plant.^9-12 MeSA is also considered a critical mobile signal for plant systemic acquired resistance (SAR).⁹ In plants, MeSA is synthesized by salicylic acid carboxyl methyltransferase (SAMT), catalyzing the transfer of the methyl group from S-adenosyl-L-methionine (SAM) to the salicylic acid (SA) carboxyl group,¹³ while MeBA is generated by SAM and benzoic acid (BA) catalyzed by BAMT. Previous studies showed that some SAMTs also had BAMT activity, and vice versa, and these enzymes were named BSMT.¹⁴ SAMTs, BAMTs and BSMTs belong to a larger family of methyltransferases called SABATH methyltransferases, which are named according to the first three enzymes (SAMT, BAMT and theobromine synthase) to be isolated and characterized in this family.¹⁵ In addition to SAMTs, BAMTs and BSMTs, the SABATH family also contains jasmonic acid carboxyl methyltransferase (JAMT/JMT), farnesoic acid MT (FAMT), indole-3-acetic acid MT (IAMT), gibberellic acid MT (GAMT), anthranilic acid MT (AAMT), cinnamate/p-coumarate MT (CCMT) and loganic acid MT (LAMT).^16-19 Additionally, except for the carboxyl MTs, the nitrogen-methyltransferases involved in caffeine biosynthesis, such as tea caffeine synthase (TCS) and theobromine synthase (TS), also belong to this family.¹⁵ The identification of SABATH members in different plants was mainly due to the significant role of their products in the plant defense response,^13,18,20-22 which indicated the special regulatory role of the SABATH family in plant defense.

The SABATH genes are widespread in plants, from algae to higher plants, where the number of family members has undergone a significant expansion.²³ As mentioned above, SABATH methyltransferases have many catalytic substrates, and a given gene in this family can also catalyze different substrates, which indicates the diversity of gene functions. For example, in Arabidopsis, there are 24 members in the SABATH gene family, AtBSMT1 and AtJMT are involved in plant defense responses, AtGAMT1 and AtGAMT2 regulate seed development and germination, and AtIAMT is involved in leaf development.^13,15 Overexpression of AtJMT gene in rice inhibited the development of rice spikelets, which led to the decrease of grain yield.²⁴ In rice, the SABATH family contains 41 members, and among the three significantly upregulated members after fall armyworm feeding, OsSABATH3 was the only member that exhibited both SAMT and BAMT activity.^3,25 Overexpression OsJMT1 in rice decreased the brown planthopper induced levels of JA and JA-Ile in transgenic plant as well as reduced plant height and yield.²² In tea plant, MeSA and MeJA not only are involved in the defense response as HIPVs but also are the components of the special aroma of tea.^26-29 Moreover, in addition to being a component of tea flavor and quality, caffeine content is also one of the indexes to evaluate the disease resistance of tea plant.¹ MeSA, MeJA and caffeine are all the products of SABATH methyltransferases, which suggests that the SABATH family in tea plant has the dual function of regulating both tea quality and tea plant defense responses. However, the gene numbers and potential functions of SABATH methyltransferases in tea plant are still unknown. Considering their special roles in the formation of HIPVs and the diversity of their gene functions, it is necessary to systematically investigate the gene numbers and underlying resistance regulatory mechanisms in tea plant. In this study, based on published genomic information of tea plant,^30,31 genome-wide identification of SABATHs was performed, and then phylogenetic analysis was used to explore the relationships of CsSABATHs with SABATHs from different plants. Finally, based on published transcriptome data and qRT-PCR, the expression patterns of CsSABATHs in different tissues and in response to MeJA, SA, JA and MeSA treatments and Ectropis obliqua damage were also analyzed. Our results provide a preliminary understanding of the potential functions of SABATHs in tea plant and will be helpful for future targeted studies.

Materials and methods

Identification of putative SABATH proteins in tea plant

To obtain putative SABATH protein sequence information for tea plant, the predicted protein sequences in the tea plant (Camellia sinensis var. sinensis) genome were first downloaded from the Tea Plant Information Archive (TPIA, http://tpia.teaplant.org).^31,32 Then, the amino acid sequences of the 24 AtSABATH proteins of Arabidopsis were obtained from The Arabidopsis Information Resource (TAIR, https://www.arabidopsis.org/) used as queries in local BLASTp searches (E-value <10⁻¹⁰) against the tea plant protein database. The obtained nonredundant candidate proteins were then confirmed in the Pfam database (http://pfam.xfam.org/) and SMART database (http://smart.embl.de/). The sequences with the PF03492 domain were considered members of the SABATH gene family of tea plant.

Sequence characteristics of the tea plant SABATH gene family

The DNA sequences and CDSs of SABATH genes in tea plant were also obtained from the TPIA. The theoretical isoelectric point (pI), molecular weights (MW) and grand average of hydropathicity (GRAVY) values of the tea plant SABATH proteins were predicted using the online tool on the ExPASy server (https://web.expasy.org/protparam/). The subcellular localization of each protein was predicted using the WoLF PSORT (https://wolfpsort.hgc.jp/) and Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/) online tools.³³

Gene structure, conserved domains and phylogenetic analysis

The exon-intron structures of each gene were analyzed using the DNA sequences and corresponding open reading frame (ORF) sequences and were visualized with the Gene Structure Display Server (GSDS) program (http://gsds.cbi.pku.edu.cn/index.php).³⁴ The conserved motifs of each protein were predicted using the Multiple Em for Motif Elicitation (MEME) online tool (http://meme-suite.org/tools/meme), and the maximum motif number was set to 20. The conserved motifs were visualized with the TBtools toolkit.³⁵ The complete amino acid sequences of CsSABATHs were aligned using ClustalX (1.83), and a phylogenetic tree was constructed using the neighbor-joining (NJ) method in MEGA6.0 with 1000 bootstrap replicates. A phylogenetic tree of the SABATH proteins from different plants was generated with the aligned sequences using the Maximum-Likelihood (ML) method in PhyML 3.0 (http://atgc.lirmm.fr/phyml/) with the Jones, Taylor and Thornton (JTT) amino acid substitution model,³⁶ and the tree was visualized using the iTOL (https://itol.embl.de/) online tool.

Stress-related cis-element identification in the promoter regions

The 2 kb upstream sequences of the start codon of each gene were considered the promoter regions. Then, the putative cis-regulatory elements were analyzed in the PlantCARE program (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).³⁷

Expression pattern analysis based on the transcriptome data

The RNA-seq data of gene expression in apical buds, young leaves, mature leaves, old leaves, stems, flowers, fruits and roots of tea plant and tea plant in response to methyl jasmonate (MeJA) treatment were downloaded from the TPIA database.³² The RNA-seq data of tea leaves in five different months (SRP119740)³⁸ and the tea plant in response to Ectropis obliqua damage (PRJNA439206)³⁹ were downloaded from the NCBI Sequence Read Archive (SRA) database. The HISAT, StringTie and Ballgown tools were used to analyze the RNA-seq data,⁴⁰ and the gene expression levels were calculated using fragments per kilobase million (FPKM) values. The accumulation data of catechins in the above eight tissues were also obtained from TPIA database, and the Gene Expression and Metabolite accumulation Correlation Analysis Tool (http://tpia.teaplant.org/Gene2Metabolite.html) was used to analyze the correlation between CsSABATHs expression and catechins content. The expression levels of CsSABATHs and the correlation relationships were visualized using the TBtools toolkit.³⁵

Plant materials, growth conditions and treatments

For the hormone treatment, one-year-old cutting seedlings of the tea plant cultivar ‘Qiancha 1ʹ (Camellia sinensis var. sinensis) were used as materials. Each seedling is about 40 cm high and contains 8–10 healthy leaves. The seedlings were collected from the nursery base of the Tea Research Institute, Guizhou Academy of Agricultural Sciences (Meitan, Zunyi, China), and were planted in pots at the greenhouse of the Tea Research Institute, Guizhou Academy of Agricultural Sciences (Guiyang, China) for hydroponics. The formulation of the hydroponic nutrient solution was as follows: (NH₄)₂SO₄, 0.713 mM; KH₂PO₄, 0.1 mM; K₂SO₄, 0.46 mM; CaCl₂, 0.5 mM; MgSO₄, 0.41 mM; Al₂(SO₄)₃ · 18H₂O, 0.1 mM; FeSO₄ · 7H₂O, 0.032 mM; EDTA-Na₂ · 2H₂O, 0.032 mM; CuSO₄ · 5H₂O, 0.002 mM; MnSO₄·H₂O, 0.09 mM; Na₂MoO₄ · 2H₂O, 0.0026 mM; ZnSO₄ · 7H₂O, 0.0091 mM; H₃BO₃, 0.046 mM; and pH4.2. The conditions of the greenhouse were as follows: 25°C/20°C (14 h light/10 h dark) and 80% relative humidity. The nutrient solution was exchanged every 3 days. After 2 weeks, the seedlings were subjected to MeSA, SA and JA treatments. For each treatment, 500 mL distilled water with 1 mM MeSA, SA or JA and 2% absolute ethyl alcohol was sprayed on the leaves of 200 tea plant seedlings. The seedlings sprayed with distilled water with 2% absolute ethyl alcohol were used as controls. One bud and two leaves from at least 10 seedlings were mixed together as one sample. The samples were collected at 0, 2, 12, 24 and 48 h after treatment. All samples were frozen in liquid nitrogen immediately and stored at −80°C for RNA extraction.

For the Ectropis obliqua feeding treatment, the two-year-old potted seedlings of ‘Qiancha 1ʹ were used as materials. The treatment of Ectropis obliqua feeding was according to the description of Yang et al. (2019).³⁹ In brief, the three instar larvae were placed evenly on the young leaves of the tea seedlings after starvation for 2 hours, and 20 insects were placed in each pot of tea seedlings. When about 1/3 of the single leaf was eaten, the insects were removed, and the samples in three biological replicates (three leaves as one biological repeat) were collected at 3 and 12 h after the feeding treatment. The seedlings without insects were used as a control. All samples were frozen in liquid nitrogen immediately and stored at −80°C for RNA extraction.

RNA isolation and quantitative real-time PCR (qRT-PCR) expression analysis

For RNA isolation, the leaves of one sample were first ground in liquid nitrogen. Then, total RNAs were extracted using the TaKaRa MiniBEST Plant RNA Extraction Kit (TaKaRa) with DNase I treatment to remove genomic DNA, and reverse transcription was performed using the PrimeScript^TM RT Reagent Kit (TaKaRa) according to the manufacturer’s instructions. The primers were designed with Primer Premier 5 and synthesized by Sangon Biotech (Shanghai, China). The 10 μL reaction contained 5 μL of GoTaq® qPCR Master Mix (Promega), 2 μL of RNase-free H₂O, 1 μL of each primer and 1 μL of cDNA. The qRT-PCR was performed with the ABI Vii ATM 7 Real-time PCR System. The CsTBP and CsActin gene were used as internal control for hormone treatment⁴¹ and insects damage treatment, respectively. Each sample was repeated three times and the relative expression levels of each gene were calculated using the 2^−ΔΔCT method.⁴² The primers used for qPCR analysis are listed in Table S1.

Results

Identification and sequence analysis of the SABATH proteins in tea plant

To fully identify the SABATH proteins in tea plant, a total of 24 SABATH family members from Arabidopsis were used as queries to perform local BLASTp (E-value <10⁻¹⁰) searches against the tea plant protein database. All the nonredundant candidate proteins were then subjected to domain analysis in the Pfam and SMART databases to confirm the presence of the PF03492 domain (Figure S1). Finally, a total of 32 SABATH proteins were identified in tea plant. These proteins were then named CsSABATH1 to CsSABATH32 according to the order of gene IDs (Table 1). According to their genomic position, most of the CsSABATH genes were distributed in different linkage groups as gene clusters (Table 1). Among them, linkage group 1 contained 6 members, linkage group 3 contained 4 members, linkage groups 5 and 6 each contained 3 members, and linkage groups 2 and 8 each contained 2 members. CsSABATH11 to CsSABATH14 and CsSABATH24 to CsSABATH27 were located on Scaffold3862 and Scaffold155, respectively. In addition, CsSABATH31, CsSABATH7, CsSABATH10 and CsSABATH9 were mapped onto linkage group 4, linkage group 9, linkage group 14 and Scaffold3313, respectively. The amino acid (AA) numbers of these proteins ranged from 142 (CsSABATH13) to 486 (CsSABATH3). The molecular weights (MW) varied from 15.5 kDa to 54.4 kDa, with an average of 35.1 kDa. The theoretical isoelectric points (pI) ranged from 4.79 (CsSABATH27) to 8.5 (CsSABATH7), with an average of 5.72. The grand average of hydropathicity (GRAVY) varied from −0.538 to 0.344, and only the GRAVY values of CsSABATH24, CsSABATH9 and CsSABATH7 were greater than 0, which indicated that the majority of the CsSABATHs were hydrophilic proteins. Subcellular localization prediction analysis showed that the majority of the CsSABATHs were located in the cytoplasm, while some of them might be in the nucleus.

Table 1. The sequence characteristics of the identified SABATH genes in tea plant.

Name	Gene ID	Map position	ORF (bp)	Animo acids (aa)	MW (kDa)	pI	GRAVY	Subcellular localization	Best match in Arabidopsis
CsSABATH1	TEA003157	LG06:Scaffold716:185070–186188	1119	372	41.5	6.89	−0.465	Cytoplasm	AT1G68040
CsSABATH2	TEA003158	LG06:Scaffold716:164662–182918	819	272	30.3	5.60	−0.538	Cytoplasm/Nucleus	AT1G68040
CsSABATH3	TEA003159	LG06:Scaffold716:192728–203176	1461	486	54.4	5.24	−0.162	Cytoplasm/Nucleus	AT1G68040
CsSABATH4	TEA004178	LG05:Scaffold5024:541830–544428	600	199	22.7	4.81	−0.151	Cytoplasm	AT1G68040
CsSABATH5	TEA007272	LG03:Scaffold2121:1751373–1754147	1041	346	38.2	5.63	−0.082	Cytoplasm	AT5G37990
CsSABATH6	TEA010054	LG03:Scaffold11064:165061–171926	555	184	21.1	4.84	−0.111	Cytoplasm	AT1G19640
CsSABATH7	TEA010123	LG09:Scaffold52:45378–60043	1011	336	37.0	8.50	0.344	Chloroplast/Nucleus	AT1G19640
CsSABATH8	TEA010250	LG05:Scaffold6545:447–23473	1011	336	37.7	5.28	−0.063	Cytoplasm/Nucleus	AT1G15125
CsSABATH9	TEA012581	Scaffold3313:1118267–1120922	1056	351	39.1	7.55	0.088	Cytoplasm	AT3G11480
CsSABATH10	TEA013347	LG14:Scaffold543:186929–195742	1275	424	47.0	5.19	−0.305	Cytoplasm	AT1G68040
CsSABATH11	TEA014697	Scaffold3862:820395–836832	846	281	31.4	5.24	−0.130	Cytoplasm	AT1G15125
CsSABATH12	TEA014699	Scaffold3862:901325–904062	1071	356	39.9	5.63	−0.220	Cytoplasm	AT1G68040
CsSABATH13	TEA014703	Scaffold3862:590152–592061	429	142	15.5	5.06	−0.065	Cytoplasm	AT5G37970
CsSABATH14	TEA014705	Scaffold3862:542784–548558	1074	357	40.0	5.37	−0.139	Cytoplasm	AT1G15125
CsSABATH15	TEA015176	LG03:Scaffold1951:451397–471786	1113	370	41.6	5.39	−0.073	Cytoplasm	AT3G11480
CsSABATH16	TEA015791	LG01:Scaffold3514:590746–596930	1110	369	41.3	5.21	−0.131	Cytoplasm/Nucleus	AT3G11480
CsSABATH17	TEA017731	LG01:Scaffold713:750784–754501	1203	400	45.5	5.21	−0.386	Cytoplasm/Nucleus	AT1G15125
CsSABATH18	TEA021081	LG03:Scaffold2268:33311–34477	786	261	29.4	5.26	−0.252	Cytoplasm	AT1G68040
CsSABATH19	TEA021329	LG01:Scaffold1968:510030–518150	1407	468	52.7	5.44	−0.215	Cytoplasm	AT1G15125
CsSABATH20	TEA021334	LG01:Scaffold1968:550286–552109	714	237	26.5	5.38	−0.042	Cytoplasm	AT1G68040
CsSABATH21	TEA022559	LG08:Scaffold502:1378140–1380475	1008	335	37.6	5.61	−0.106	Cytoplasm/Nucleus	AT3G11480
CsSABATH22	TEA022575	LG08:Scaffold502:1263320–1265532	735	244	27.3	5.16	−0.009	Cytoplasm/Nucleus	AT3G11480
CsSABATH23	TEA024636	LG02:Scaffold656:435703–436328	459	152	16.9	5.71	−0.179	Cytoplasm	AT3G11480
CsSABATH24	TEA028049	Scaffold155:534895–536130	660	219	24.4	8.10	0.064	Cytoplasm/Nucleus	AT4G36470
CsSABATH25	TEA028050	Scaffold155:356664–359057	1005	334	37.4	5.44	−0.045	Cytoplasm/Nucleus	AT3G11480
CsSABATH26	TEA028051	Scaffold155:450150–452541	1005	334	37.8	6.02	−0.070	Cytoplasm/Nucleus	AT3G11480
CsSABATH27	TEA028052	Scaffold155:537090–543289	480	159	18.0	4.79	−0.356	Cytoplasm	AT1G68040
CsSABATH28	TEA030024	LG02:Scaffold1313:407963–410010	663	220	24.9	5.40	−0.187	Cytoplasm/Nucleus	AT3G11480
CsSABATH29	TEA031418	LG01:Scaffold3918:1284776–1292837	1077	358	40.1	5.31	−0.277	Cytoplasm/Nucleus	AT1G19640
CsSABATH30	TEA031962	LG05:Scaffold982:1605217–1609749	1164	387	42.4	5.42	−0.104	Cytoplasm/Nucleus	AT5G55250
CsSABATH31	TEA032424	LG04:Scaffold2907:173616–179621	1176	391	43.0	5.81	−0.126	Cytoplasm/Nucleus	AT5G55250
CsSABATH32	TEA032708	LG01:Scaffold1618:1069456–1075932	1077	358	40.9	7.57	−0.180	Cytoplasm/Nucleus	AT5G66430

Gene structure, conserved domains and phylogenetic analysis of the CsSABATH gene family

To investigate the structural evolution of CsSABATH genes in tea plant, the intron-exon structure of each gene was analyzed based on the DNA sequence and CDS. In addition, an NJ evolutionary tree was constructed to explore the relationships between exon-intron distribution patterns and phylogenetic classification. Our results showed that the 32 CsSABATHs could be divided into three groups (I, II and III), the group I and group III contained the same number (15) of members, while group II only contained two members (CsSABATH30 and CsSABATH31) (Figure 1(a)). The intron numbers of CsSABATH genes varied from 0 (CsSABATH1) to 6 (CsSABATH3), and the majority of them (27/32) had no more than three introns (Figure 1(b)). The CsSABATHs clustered onto the same branch and always had similar intron-exon structures. Then, the conserved motifs of each protein were further searched in MEME. As shown in Figure 1(c), the CsSABATH proteins in the same subgroups also had similar motif types and distribution patterns. Motif 1 to 7 existed in almost all the CsSABATH proteins, while motif 8 and 10, motif 19 and 9, and motif 12 and 13 were specifically distributed in the group I, II and III, respectively.

Figure 1. Phylogenetic, gene structure and conserved motif analyses of CsSABATHs. (a) Phylogenetic relationships of the 32 CsSABATH proteins. (b) Exon-intron structures of CsSABATH genes. (c) The conserved motifs of CsSABATH proteins.

Phylogenetic analysis of the SABATH gene family in plants

To further clarification the evolutionary relationship and to predict the potential functions of CsSABATH genes, the genome-wide SABATH proteins from Marchantia polymorpha, water lily, kiwifruit, cacao, coffea, grape, Arabidopsis, rice, Picea abies and Camellia sinensis var. assamica were identified and downloaded from public databases (Table S1). Additionally, the SABATH proteins whose functions have been reported were also added. Finally, a total of 308 SABATH proteins from 19 plant species were used to reconstruct the phylogenetic tree. As shown in Figure 2, these proteins could also be divided into three major groups (group I, group II and group III), and group I contained the most SABATH members. The N-methyltransferases of cacao, coffea and tea plant were mainly clustered in group I but in three distinct clades. In addition, carboxyl-methyltransferases such as AAMT, BAMT, SAMT, BSMT and JMT were also clustered in group I. Among them, CsSABATH32 was clustered together with the SAMT and BSMT genes of other plants, while CsSABATH29 had a close relationship with AtJMT and OsJMT1. The SABATH members of Camellia sinensis var. assamica were also mainly concentrated in group I. Although group II included the smallest number of SABATH members, it contained the most plant species, and the 12 SABATH gene family members of Marchantia polymorpha were all clustered in this group. Additionally, the carboxyl-methyltransferases JMT, GAMT, CCMT and IAMT were also clustered in this group, and CsSABATH30 and CsSABATH31 had a close relationship with IAMTs. The members of group III included only dicotyledon SABATH genes. With the exception of FAMTs, most of the members in group III were functionally unknown. A total of 15 SABATH genes of Camellia sinensis var. sinensis (CSS) were clustered in group III, while only two SABATH genes of Camellia sinensis var. assamica (CSA) were clustered in this group, and the group III SABATH genes of tea plant displayed clear expansion in the CSS genome (Figure 2, Figure S2).

Figure 2. Phylogenetic analysis of SABATH proteins in different plants. The accession number of each protein is listed in Table S2.

Analysis of cis-elements in the promoter regions of CsSABATHs

Cis-elements in the promoter region always play a regulatory role in gene functions. Here, the 2000 bp upstream sequence from the start codon of each CsSABATH gene was obtained, and the cis-elements were then predicted in PlantCARE. The results showed that all the candidate promoter regions possessed the typical TATA boxes and CAAT boxes. Then, we analyzed the regulatory elements associated with stress and hormonal responses. As shown in Figure 3, a total of 14 kinds of stress or hormonal response cis-regulatory elements were identified, including anaerobic induction element ARE, defense and stress responsiveness element TC-rich repeats, low-temperature responsiveness element LTR, drought inducibility element MBS, ABA responsiveness element ABRE, auxin-responsive element AuxRR-core and TGA-element, gibberellin responsive element P-box, GARE-motif and TATC-box, MeJA responsiveness element CGTCA-motif/TGACG-motif and SA responsiveness element TCA-element. The element types and numbers of each gene varied, which indicated their functional differences.

Figure 3. The statistics of the major stress-related and hormonal response cis-regulatory elements detected in the promoter regions of CsSABATH genes.

Expression patterns of the CsSABATH genes revealed by transcriptome analysis

Based on public transcriptome data, the expression of CsSABATH genes in different tissues of tea plant and different seasons of tea leaves was used to analyze their spatiotemporal expression patterns. The expression of 27 CsSABATH genes was detected in at least one of the eight tea plant tissues, and the expression levels showed diversity. As shown in Figure 4(a), the majority of the CsSABATH genes showed tissue-specific expression patterns. For example, CsSABATH9, 15, 16, 24, 25, 28 and 29 were mainly expressed in young leaves and young stems, CsSABATH26 and 32 showed high expression in flowers, and the other 18 members except CsSABATH7 were highly expressed in roots. CsSABATH7 showed very little expression in all eight tissues but showed specific expression in the stem. Among the 27 CsSABATH genes, CsSABATH16 and CsSABATH25 were found to have extremely high expression in all eight tissues compared with that of the other members. In different seasons, the CsSABATH genes also exhibited a variety expression patterns (Figure 4(b)). Among them, CsSABATH4, 9, 12, 15, 16, 17, 18, 23, 24, 25, 26 and 28 showed high expression in the leaves of April, CsSABATH1, 2, 3, 5, 7, 10, 19 and 20 were mainly expressed in June, while CsSABATH8, 14 and 31 were specifically expressed in tea leaves only in October.

Figure 4. The spatiotemporal expression patterns of SABATH genes in tea plant. (a) The expression patterns of CsSABATHs in eight different tissues of tea plant. (b) The expression patterns of CsSABATHs in different seasons of tea leaves.

Since the promoter analysis showed that many cis-elements associated with low-temperature, drought, MeJA and defense and stress responsiveness, we examined the expression of CsSABATHs under cold acclimation, PEG and NaCl treatments, MeJA treatment and Ectropis obliqua damage based on RNA-seq data. However, compared with MeJA and Ectropis obliqua feeding treatments, CsSABATH genes showed no obvious response to low temperature and drought or salt stress (Figure 5), which indicated that they might not have special functions in response to abiotic stress. Then, we focused on their response to the latter two treatments. In the MeJA treatment, a total of 30 CsSABATH genes were detected the expression (Figure 5(d)). Among them, CsSABATH7, 11, 12, 16, 21, 25, 26, 28, 29 and 32 were upregulated in expression after MeJA treatment, while the other 20 members showed downregulated expression. Among the members with upregulated expression, CsSABATH7, 25, 26 and 29 showed peak expression 24 h after treatment, and the other 6 members reached the peak at 12 h, then gradually decreased, and most of them exhibited recovery to the pretreatment level at 48 h. Among the members with downregulated expression, CsSABATH9, 17, 27 and 31 showed slightly upregulated expression 12 h after treatment, while CsSABATH24 and 30 began showing upregulated expression at 48 h. In the Ectropis obliqua damage treatment, a total of 25 CsSABATH genes were detected the expression (Figure 5(e)). Among them, a total of 18 members (CsSABATH1, 3, 5, 7, 8, 11, 12, 14, 15, 16, 18, 24, 25, 26, 28 and 32) exhibited upregulated expression after Ectropis obliqua fed on tea leaves, especially for CsSABATH11, 14, 16, 24, 25 and 32, which indicated the special role of them in response to the feeding of Ectropis obliqua.

Figure 5. The stress response expression patterns of CsSABATHs in tea leaves. (a) The expression patterns of CsSABATHs in response to cold acclimation. The five stages represent nonacclimated at 25 ~ 20°C (CK), fully acclimated at 10°C for 6 h (CA1-6 h) and 10 ~ 4°C for 7 days (CA1-7d), cold response at 4 ~ 0°C for 7 days (CA2-7d) and recovering under 25 ~ 20°C for 7 days (DA-7d). (b) The expression patterns of CsSABATHs in response to 25% PEG treatment for 0, 24, 48 and 72 h. (c) The expression patterns of CsSABATHs in response to 200 mM NaCl treatment for 0, 24, 48 and 72 h. (d) The expression patterns of CsSABATHs in response to 0.25% (v/v) water solution of MeJA treatment for 0, 12, 24 and 48 h. (e) The expression patterns of CsSABATHs after Ectropis obliqua feeding in tea leaves.

Based on the expression levels of different genes in the eight tissues of tea plant mentioned above, the expression correlations between CsSABATH genes and MYB, ARF, NAC, SPL, EIL, AP2/ERF and WRKY transcription factors were analyzed to explore the transcriptional regulation of CsSABATH genes. The coexpression network of them is shown in Figure 6. The expression of many transcription factor genes and CsSABATH gene family members showed significant correlations in different tissues of tea plant (Table S2), which indicated that the expression of CsSABATH genes in tea plant might also be regulated by different transcription factor genes.

Figure 6. Transcriptional regulation of CsSABATH genes. The coexpression network of CsSABATH genes with MYB, ARF, NAC, SPL, EIL, AP2/ERF and WRKY transcription factors was constructed based on their expression patterns in eight different tissues of tea plant in Cytoscape v3.7.1. The expression correlations between CsSABATH genes and different transcription factors are shown with gray lines (Pearson’s correlation test, P ≤ 0.05), and their detailed correlations are listed in Table S3.

Meanwhile, based on the expression levels of CsSABATH genes in different tissues of tea plant and the accumulation of different catechins in the tissues, the correlation of them was further analyzed. As shown in Figure S3, the majority of CsSABATH genes showed a positive correlation with the nonsoluble proanthocyanidins (PAs), while CsSABATH29, 25, 15, 24, 16 and 28 showed a strong correlation with the accumulation of total catechins, galloylated flavan-3-ols, cis-flavan-3-ols, epicatechin-3-gallate (ECG), epigallocatechin-3-gallate (EGCG) and soluble PAs.

Expression patterns of the CsSABATH genes in response to JA, SA, MeSA and Ectropis obliqua feeding

To further explore the functions of CsSABATH genes, a total of eight members from different groups based on phylogenetic relationships were selected to analyze their expression patterns under JA, SA, MeSA treatments and Ectropis obliqua feeding by qRT-PCR. As shown in Figure 7, the eight genes showed differential response patterns to the different phytohormone treatments. To our surprise, almost all of them (except CsSABATH13) were obviously downregulated in response to JA and SA treatments, and the downregulation trend lasted at least 48 h. Under JA treatment, CsSABATH13 expression was repressed at 12 h, similar to the pattern observed for other members, but then was induced peaked at 24 h, and was downregulated again. Under MeSA treatment, CsSABATH16, CsSABATH29 and CsSABATH31 were significantly upregulated at 2 h, while CsSABATH27 was significantly downregulated at this time point. However, 12 h after treatment, all of them except CsSABATH23 were downregulated in expression, and the downregulation trend of CsSABATH16 continued until 48 h after treatment. Similar to the result in the JA treatment, CsSABATH13 expression was also significantly upregulated 24 h after MeSA treatment.

Figure 7. The relative expression levels of CsSABATHs at different time points after JA, SA and MeSA treatments in tea leaves. The bars represent the mean+SE (n = 3). Significant differences between the treatments and CK were determined by Student’s t test (*p < .05, **p < .01).

The expression patterns of these CsSABATH genes in response to Ectropis obliqua feeding were shown in Figure 8. We didn’t detect the expression of CsSABATH27. Among the remaining seven members, CsSABATH8, 13, 16, 29 were significantly upregulated, while CsSABATH31 was significantly downregulated at 12 h after the damage of Ectropis obliqua.

Figure 8. The relative expression levels of CsSABATHs at different time points after Ectropis obliqua feeding in tea leaves. The bars represent the mean+SE (n = 3). Significant differences between the treatments and CK were determined by Student’s t test (*p < .05, **p < .01).

Discussion

SABATH methyltransferases are a multifunctional gene family that are ubiquitous in higher plants, and both their substrates and products play an important role in many biological processes in plants.^13,23,43-45 In this study, a total of 32 SABATH genes were identified in the tea plant genome, which was larger than the numbers in Arabidopsis (24),¹⁵ Populus (28),^46,47 Salvia miltiorrhiza (30),⁴⁸ Picea abies (10)⁴⁹ and water lily (13),⁵⁰ but smaller than the number in rice (41).³ Genome location annotation revealed that most CsSABATHs were distributed on different linkage groups of the tea plant genome in the form of gene clusters, and phylogenetic analysis showed that the members of the same gene cluster gathered on the same evolutionary branch, which indicated that the members of the same gene cluster came from a common ancestor. This phenomenon has been observed in Populus,⁴⁷ suggesting that tandem duplication is also responsible for the expansion of SABATH genes in tea plant. However, the members from the same gene clusters exhibited different expression patterns. For example, in the CsSABATH1-CsSABATH3 gene cluster, the tissue expression pattern analysis showed that CsSABATH1 was especially expressed in the root of tea plant, while CsSABATH3 was mainly expressed in the root and stem, with a higher expression level in the stem than in the root. Under MeJA treatment, the expression of CsSABATH1 was downregulated, while that of CsSABATH3 was upregulated. Although the expression of both of them was induced after Ectropis obliqua feeding on tea leaves, the expression of CsSABATH3 was upregulated to a much greater degree than that of CsSABATH1. In the above different tissues and treatments, the expression of CsSABATH2 was not detected, indicating that the function of CsSABATH2 may have been lost. These results suggested that functional divergence occurred among the three genes in this gene cluster. Additionally, functional divergence among members of the same gene cluster was also observed for other clusters of CsSABATH genes, which was indicative of the complexity of this family’s functions in tea plant.

In this study, both the NJ method and ML method were used to construct phylogenetic trees of CsSABATH proteins, and the 32 CsSABATH proteins were divided into three groups. The three class designations were also supported by gene structures and conserved motif analysis. Additionally, except for the SABATH genes of tea plant, the SABATH family members of other higher plants also contained three distinct clades. The group II SABATH family in tea plant contained only two members, CsSABATH30 and CsSABATH31, which were clustered together with IAMTs from rice, Arabidopsis and Poplar. Previous studies showed that compared with other SABATH family members showing multisubstrate catalytic activity, AtIAMT1, OsIAMT1 and PtIAMT1 all showed catalytic activity specific to IAA.^3,51 Moreover, the conserved sequence and function of IAMT in monocotyledonous and dicotyledonous plants indicated that IAMT appeared before the differentiation of dicotyledonous and monocotyledonous plants, and IAMT is thus probably an older member of the SABATH family.³ In this study, a total of 12 SABATH members were identified in the Marchantia polymorpha genome, and all of them were clustered in group II, which indicated that the SABATH gene family in Marchantia polymorpha has not yet shown clear functional differentiation. As this species is a bryophyte, this result partly supported the inference that the members of group II came from the early evolution of the SABATH family. Substrate preference was the driving factor for the functional evolution of the SABATH family.⁵² Here, we found that the members of the group II also included JAMT, GAMT, SAMT, CCMT, and others, which indicated that SABATHs with similar kinships also had diverse functions in different plant species.^44,47 While we did not obtain clear phylogenetic relationships for SABATHs in higher plants based on the phylogenetic tree, we did find several interesting features. First, group I contained the most SABATH members and the majority of them were N-methyltransferases. The N-methyltransferases of cacao, coffea and tea plant in group I were on three separate branches, indicating that they had undergone independent evolution.^30,31 In addition, this group also contained AAMT, JMT, SAMT and BSMT, which was related to the fact that their ancestor (xanthine methyltransferase, XMT) had both O- and N-methylation activities.⁵³ Second, group III contained only SABATH members from dicotyledons, where the number of SABATH family members in CSS was significantly expanded compared with CSA. Previous studies showed that group III SABATH proteins had high activity toward FA, which could convert FA to MeFA.^14,47 MeFA is an analog of the insect juvenile hormone III (JHIII), and the high content of JHIII has a special role in plant defense against insects by interfering in insect development.⁵⁴ In Arabidopsis, the FAMT could induce plant defense responses at the gene expression level.¹⁴ In this study, the group III members CsSABATH8, 11 and 13 showed significant induced after feeding by Ectropis obliqua, suggesting that the group III members may play a special role in plant defense against insect herbivory. Additionally, CSS and CSA are the two main types of cultivated tea plant varieties, and they are distinct in terms of both biological characteristics and tea quality.³¹ The biological changes brought about by the obvious expansion of group III SABATH genes in CSS are worth further exploration.

In this study, we focused on the SABATH genes in tea plant, mainly because they might have the dual functions of regulating the defense response of tea plant and affecting the quality of tea. JA, SA and their methyl esters MeJA and MeSA play a special role in plant biotic defense,^11,55-59 and the methylation of JA and SA is catalyzed by JMT and SAMT, respectively. The JMT and SAMT genes in different plants have been shown to be involved in plant defense.⁶⁰ For example, in Arabidopsis, overexpression of AtJMT1 increased the release of MeJA and conferred a higher level of resistance to Botrytis cinerea compared with that in the wild type.⁶¹ In soybean, overexpression of GmSAMT1 conferred resistance in transgenic plants to soybean cyst nematode.⁶² Here, based on RNA-seq data, the expression patterns of CsSABATHs in response to Ectropis obliqua damage were analyzed. We found that the majority (18/25) of the expressed CsSABATHs were upregulated in tea leaves after the feeding of Ectropis obliqua, which indicated that SABATH genes responded to herbivore induction in tea plants. This hypothesis was further confirmed by qPCR analysis. Among the eight selected members used for qPCR, except CsSABATH27 was not detected, four of the other seven genes were significantly upregulated in the leaves of tea plants fed by Ectropis obliqua. Moreover, CsSABATH29 was not detected in transcriptome data, but it was obviously upregulated in our qPCR analysis, which may be related to the technical differences between the two assays or the different samples used. Therefore, in addition to the obvious upregulated members identified by RNA-seq (CsSABATH11,14,16, 24, 25 and 32), the new upregulated members identified by our qPCR analysis (CsSABATH8, 13 and 29) can be used as candidate genes with a role in the defense response of tea plant, which warrants further exploration. Additionally, the expression patterns of CsSABATHs in MeJA treated tea plant leaves were also examined. We found that the genes that were highly expressed in the buds and young leaves of tea plant were also obviously induced after 12 h of MeJA treatment, which indicated that the function of the gene was closely related to the specificity of its tissue expression. Meanwhile, the promoter regions of most CsSABATH genes had different numbers of MeJA responsiveness elements. However, there was no significant correlation between expression trend and the number of relevant cis-regulatory elements in the promoter regions. For example, CsSABATH29 was mainly expressed in the buds and young leaves, was significantly induced after MeJA treatment and peaked at 24 h, but its promoter region did not contain MeJA responsiveness elements. These phenomena were also observed in our JA, SA and MeSA treated gene expression analysis, which indicated that the gene functions could be influenced by cis-regulatory elements, but do not depend on these cis-regulatory elements.⁶³ Generally, MeJA and JA always exhibit similar signaling functions in plants. However, our results showed that the response of CsSABATHs under MeJA and JA treatment was not consistent. In addition to the difference between RNA-seq and qPCR detection methods, it may be related to the difference in the response mechanism of plants to MeJA and JA.⁶⁴ Previous studies have found that the resistance of plants to herbivore caused by exogenous MeJA is not directly from MeJA itself, but from its demethylated product JA.^22,65 Exogenous MeJA treatment promotes the production of endogenous JA,⁶⁵ while the member of SABATH family JMT can catalyze JA to produce MeJA, which may be one of the reasons for the difference of CsSABATHs expression between MeJA and JA treatment. In addition, phylogenetic analysis showed that CsSABATH29 had a close relationship with AtJMT and OsJMT1. Therefore, CsSABATH29 could act as a JMT in tea plant, and its special functions deserve further study.

In this study, regardless of the public transcriptome data or qPCR analysis, the expression patterns of CsSABATHs were diverse, which also indicated the complexity of SABATH genes in regulating various biological functions of tea plant. The biological processes of plants are not regulated by a single gene, but by a series of genes interacting with each other through a gene network.⁶⁶ HIPVs, such as MeSA and MeJA are all secondary metabolites. The biosynthesis of these secondary metabolites is not only catalyzed by related enzymes but also regulated by many transcription factor genes. A classic example is the MYB-bHLH-WD40 (MBW) ternary transcription factor complex, which regulates the synthesis of anthocyanins and proanthocyanidins in plants.^67,68 However, whether the synthesis of these volatile methyl esters is regulated by transcription factor genes is largely unknown. Generally, the expression patterns of synthases and their regulated transcription factor genes are significantly correlated, and many regulatory relationships have been verified based on these in different plants.^69,70 Here, transcription factor families, including MYB, bHLH, AP2/ERF, WRKY, SPL, EIL, ARF and NAC, were identified from the tea plant genome. Based on their expression patterns in eight different tissues of tea plant, Pearson’s correlation analysis of them with CsSABATH genes was performed. Our results showed that the expression patterns of many transcription factor genes were significantly correlated with different members of the CsSABATH gene family, suggesting a potential regulatory relationship between them. Moreover, we also found a clear correlation between different CsSABATH family members, and a CsSABATH gene might be regulated by multiple transcription factors, which further demonstrated the complexity of the regulation of plant secondary metabolism.

Conclusions

In this study, we identified 32 SABATH methyltransferase genes in the tea plant genome, which were divided into three groups based on intron-exon structures, conserved motifs and phylogenetic analysis. Tandem duplication might be responsible for the expansion of SABATH genes in tea plant, and functional divergence among these gens during their evolution was also observed according to expression patterns analysis. Compared with those in CSA, the group III SABATH family members in CSS were significantly expanded, but the difference in biological function was unclear. Ectropis obliqua feeding induced upregulated expression in most SABATH genes in tea plant, and CsSABATHs also presented various expression patterns in the treatments with the plant defense response factors SA, JA, MeJA and MeSA, which indicated the special roles of CsSABATHs in tea plant defense responses. The transcriptional level of CsSABATH8, 11, 16, 25, 29 and 32 was significantly induced after MeJA treatment and feeding on Ectropis obliqua, and their roles in tea plant defense response should be further explored. Additionally, the biological function of SABATH genes might be regulated by many transcription factor genes, such as MYB, bHLH, AP2/ERF, WRKY, SPL, EIL, ARF and NAC genes. The findings of this study provide a preliminary understanding of the potential functions of SABATHs in tea plant and will be helpful in their targeted functional characterization in the future.

Author contributions

Yan Guo: Methodology, Data analysis, Writing-original draft preparation. Dahe Qiao: Methodology, Visualization, Writing-review & editing, Funding acquisition. Chun Yang: Data analysis, Writing-original draft preparation. Juan Chen: Funding acquisition, Writing-review & editing. Yan Li: Data analysis. Sihui Liang: Data analysis. Kaiqin Lin: Data analysis. Zhengwu Chen: Conceptualization, Supervision, Writing-review & editing.

Disclosure of potential conflicts of interest

The authors declare that they have no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

This research was funded by the Science and Technology Foundation of Guizhou Province (grant number 20191303 and 20201Y119), the Youth Foundation of Guizhou Academy of Agricultural Sciences (grant number 201845 and 201816), the Science and Technology Innovation Platform Construction of Guizhou Province (grant number 20202103) and the Science and Technology Support Plan of Guizhou Province (grant number 20172557).