Genome-wide identification of the OMT gene family in Cucumis melo L. and expression analysis under abiotic and biotic stress

Plant and material

‘Yangjiaomi’ is widely cultivated in China and has become an important commercial cultivar because of its crunchy texture, juiciness, and sweet taste. ‘Yangjiaomi’ was growing in the plant culture room of Liaocheng University, Liaocheng, China. Melon seeds were germinated in an incubator at 28 °C and then sown in 50-hole cavity trays (grass: vermiculite: perlite = 1:1:1, v/v). After the seedlings had grown three leaves and heart-shaped leaves, the treatments were applied in a plant culture room (26 °C/20 °C, day/night) with a 16 h photoperiod. The same melon plants were selected for the abiotic stress treatments, which were set according to a previous study: NaCl stress (300 mM for 1, 3, 5 and 7 d), chilling stress (6 °C for 6 and 12 h), and high temperature and humidity (HTH) stress (day/night temperature at 45 °C/35 °C, soil moisture at 100%, and humidity at 90%) (Wang et al., 2016; Diao et al., 2020; Weng et al., 2022). The control and treated melons were collected to determine the relative expression of the CmOMT gene family. Each stress treatment group included 15 melon plants, and three biological replicates were performed for all treatment groups.

Identification of O-methyltransferase (OMT) genes in Cucumis melo L.

The melon (DHL92) protein database (Melon protein v4.0) was downloaded from the Cucurbitaceae Genome Database (CuGenDB; http://cucurbitgenomics.org/). To comprehensively identify the melon OMT gene family members, we searched the NCBI (Conserved Domain Search, http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi/) and PFAM websites (http://pfam.xfam.org/) and analyzed the conserved motifs and protein functions of the reported ASMT and COMT proteins. PF00891 was employed as the reference sequence for the Hidden Markov Model (HMM) analysis of the melon protein database, and proteins with an e-value ≤ 1 × e ⁻¹⁰ were selected as candidate OMT proteins (Liu et al., 2015; Chang et al., 2021). The candidate proteins were subjected to protein sequence identification using the SMART website (http://smart.embl-heidelberg.de/) (Letunic, Khedkar & Bork, 2017). Sequences without OMT structural domains or incomplete structural domains were removed, and the longest sequence was retained if redundant sequences were identified.

Physicochemical properties and chromosomal mapping of CmOMT genes in Cucumis melo L.

The molecular weight, isoelectric point (PI), formula, instability index, aliphatic index, and grand average hydropathicity (GRAVY) of the CmOMT genes were analyzed using the online tool ExPASy (http://web.expasy.org/protparam/) with the default parameters. The melon genome website (http://cucurbitgenomics.org/organism/18) provides the position information of CmOMT genes using MapChart software to map the distribution of CmOMT genes on chromosomes. The identification criteria for tandem duplication genes were as follows: (i) adjacent homologous genes were located on the same chromosome, with only one intercalated gene in the middle, and (ii) the comparison length and similarity of the two gene sequences were more than 70% (Zhu et al., 2014; Zhao et al., 2018).

Phylogenetic analysis, exon-intron structure, and motif analysis

Multiple alignments of the amino acid sequences of the 21 CmOMT proteins were performed using the ClustalW tool with default parameters (Larkin et al., 2007). To investigate the phylogenetic relationships of OMT homolog protein sequences among six species (Arabidopsis thaliana, Cucumis melo L., Cucumis sativus L., Solanum lycopersicum L., Oryza sativa L., and Citrullus lanatus (Thunb.) Matsum. et Nakai), a phylogenetic tree was constructed using the neighbor-joining algorithm of MEGA X software. A phylogenetic tree was constructed using the neighbor-joining Jones Taylor–Thornton (JTT) matrix model in MEGA X software, where bootstrap repetitions were set to 1,000 to assess the plausibility of the evolutionary tree (Kumar, Stecher & Tamura, 2016). The phylogenetic tree was visualized using Figtree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

The CmOMT gene structure was analyzed using the Gene Structure Display Server (GSDS) website (Hu et al., 2015) (http://gsds.cbi.pku.edu.cn). General Feature Format (GFF) annotations for melons were downloaded from the genome database. The required annotation content for the CmOMTs exon-intron structure was uploaded to the GSDS website. The online tool Multiple Em for Motif Elicitation (MEME: https://meme.nbcr.net/) was used to further identify and analyze the conserved protein motifs of the CmOMT gene family members, and a motif number of 5 was implemented (Bailey et al., 2009).

Promoter analysis, transmembrane helix prediction, and structure prediction for CmOMT proteins

The promoter sequences of the CmOMT genes were analyzed using PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). Briefly, we used the sequence 2,000 bp upstream of the ATG start codon of the CmOMT gene family as the promoter sequence. Cis-regulatory elements in CmOMT gene promoters identified using PlantCARE were conserved for visualization using the GSDS website. In standard mode, the three-dimensional (3-D) structure of each CmOMT protein was modeled using the PHYRE2 tool (http://www.sbg.bio.ic.ac.uk/phyre2/html/page.cgi?id=index).

GO term and expression analysis for CmOMT genes

Gene Ontology (GO) term analyses were performed for CmOMT genes using the R package (ClusterProfiler package, enrichplot package), and terms with Q-values ≤ 0.05 were considered significantly enriched.

The expression of CmOMT genes in various organs was examined using a previously published RNA-seq dataset (Yano et al., 2020). A heatmap was generated to illustrate the spatiotemporal expression in the callus, dry seeds, root, stem (downside and upside), shoot apex, leaves (young and 6th–12th), tendril, flower (anther male, petal female, and stigma female), ovary [0–4 DAF (days after flowering)], fruit flesh (8–50 DAF), and fruit epicarp (8–50 DAF). Published RNA-seq datasets analyzed by Sebastiani et al. (2017) and Zhu et al. (2018) were analyzed for Fusarium oxysporum f. sp. melonis race 1.2 (FOM 1.2) and powdery mildew, respectively. Published RNA-seq data were used to generate heatmaps and volcano plots using the R package (Heatmaply package) and online tool Dynamic volcanogram (https://www.omicshare.com/tools/), respectively. To verify the expression of CmOMT genes under different stress conditions, total RNA was extracted using the RNA isolator Total Extraction Reagent (Vazyme, Nanjing, China). The extracted RNA was reverse-transcribed to cDNA using the HiScript II First Strand cDNA Synthesis Kit (Vazyme, Nanjing, China). qRT-PCR assays were performed on a qTOWER³G Real-time System (Analytik Jena AG, Jena, Germany). The relative expression based on three biological and technical repeats was calculated by the 2^−ΔΔCT method (Livak & Schmittgen, 2002). All primers used for qRT-PCR are listed in Table S1.

Statistical analysis

All data were analyzed using SPSS software version 26.0 (SPSS Inc., Chicago, IL, USA) and presented as the mean ± standard deviation (SD) of three biological repeats. Differences among the results were tested by one-way ANOVA (one-way analysis of variance), with a P-value ≤ 0.05 indicating significance.

Identification and characterization of CmOMT genes in Cucumis melo L.

Twenty-one CmOMT candidate genes were identified in the melon genome (Table S2) and named CmOMT1 to CmOMT21 based on the chromosomal gene order. Certain CmOMT genes were also distributed in adjacent regions on the same chromosome (Table 1 and Fig. 1A). The genes presented an uneven distribution on 12 chromosomes and were present on six chromosomes; however, four CmOMTs were not mapped to any of the chromosomes (Fig. 1A). CmOMTs were intensively distributed on chromosome 1, three CmOMTs were distributed on chromosome 6, two CmOMTs were distributed on chromosomes 4 and 11, and only one CmOMT gene was distributed on chromosomes 7 and 10 (Fig. 1A). According to a previously defined tandem arrangement, CmOMT1–CmOMT4, CmOMT13, CmOMT14, and CmOMT18–CmOMT21 were not tandem duplicated genes, whereas other CmOMTs appear to be generated by tandem repeat events (Zhu et al., 2014; Zhao et al., 2018). For example, CmOMT5–CmOMT8 may belong to a tandem duplication gene cluster while CmOMT9–CmOMT12 may be part of another tandem duplication gene cluster.

The chromosome location, strand, description, number of amino acids, molecular weight, PI, formula, instability index, aliphatic index, and grand average hydropathicity (GRAVY) characteristics of the CmOMT genes were obtained. The average amino acid sequence of the 21 CmOMT proteins was 298 amino acids, with the value ranging from 93 (CmOMT4) to 384 amino acids (CmOMT8) (Table 1). The average molecular weight of the 21 CmOMT proteins was 33,211.24, the PI was approximately 5.50, the instability index was approximately 38.88, and the aliphatic index was approximately 91.63 (Table 1). In addition, five of the CmOMT proteins (CmOMT3, CmOMT4, CmOMT9, CmOMT15, and CmOMT18) were hydrophobic (Table 1).

Analysis of the gene structure and motifs in CmOMT genes

To further investigate the CmOMT gene family, the UTR and CDS distribution and number of CmOMT genes were analyzed based on phylogenetic trees. On the same branch, the CmOMTs were divided into class I and class II, which contained nine and six CmOMTs, respectively (Fig. 1B). On the other branch, six CmOMTs were grouped into class III, which was distinct from classes I and II (Fig. 1B). Genetic structural characterization showed that the number of CDS in classes I and II was similar, with 1–3 CDS regions (Fig. 1B). Class III differed from classes I and II in terms of gene structure and included 4 CDS regions, excluding CmOMT18 (Fig. 1B). Each exon in the same sister gene pair was very similar and presented the same size. However, introns were more variable in sequence length than the strongly conserved exons of sister gene pairs. This result indicates that introns are evolutionarily more unstable than exons, which is consistent with traditional evolutionary theory (Zhang et al., 2021a; Zhang et al., 2021b).

Conserved protein motifs were identified in CmOMTs. Five different motifs with 15–34 amino acids were found in the CmOMT gene family (details on the conserved structural domains are illustrated in Fig. S1A). CmOMT protein sequences were compared using a phylogenetic tree, and the locations were mapped using clustering and motif analyses. Details of the CmOMTs are presented in Fig. S1B. Motif 1 was present in most of the CmOMTs except for CmOMT2 (Fig. 1B). Motif 2 was not present in CmOMT3, CmOMT4, CmOMT7, CmOMT12, and CmOMT15 (Fig. 1B). Motifs 3 and 4 represented the basis of the CmOMT structural domain because all CmOMTs contained motifs 3 and 4. Motif 5 was identified in only nine CmOMT members, which belonged to class I (CmOMT7, CmOMT8, CmOMT13, and CmOMT19) and class II (CmOMT9–CmOMT12 and CmOMT21) (Fig. 1B). These results suggested that there were significant differences between CmOMT proteins after whole genome duplication (WGD) events and that the loss and retention of motifs may be related to the evolution of subbranches.

Phylogenetic tree and structure analysis of OMT proteins

To investigate the evolutionary relationships of plant OMT homologs, 120 OMTs were identified in melon, Arabidopsis (Chen et al., 2021), cucumbers (Table S3), tomatoes (Lu et al., 2019), rice (Liang et al., 2022) and watermelon (Chang et al., 2021), and an unrooted phylogenetic tree was constructed using the NJ algorithm. The phylogenetic tree branches showed that OMTs from the selected plant species could be divided into class I, class II, and class III, which contained 53, 44, and 23 OMTs, respectively (Fig. 2). Class I contained nine CmOMTs, six CsOMTs, ten SlOMTs, 24 OsOMTs, and four ClOMTs, and it was not included among the AtOMTs (Fig. 2). Class II contained six CmOMTs had 15 AtOMTs, six CsOMTs, four SlOMT genes, nine OsOMTs, and four ClOMTs (Fig. 2). Class III did not contain OsOMTs and was composed of six CmOMTs, two AtOMTs, five CsOMTs, two SlOMTs, and eight ClOMTs (Fig. 2). Multiple species form gene clusters at the base of the NJ tree, with melon, cucumber, tomato, and watermelon showing a closer homology of OMT proteins than Arabidopsis and rice, which may reflect the functional diversification of the OMT gene family after the evolution of dicotyledonous and monocotyledonous plants. In addition, cucurbit OMT proteins (CmOMTs, CsOMTs, and ClOMTs) showed a close homology, suggesting that the OMT family experienced species-specific amplification during the evolutionary process.

The structure of a protein determines its function. Therefore, the 3-dimensional structures of CmOMTs were predicted using the PHYRE2 tool using the standard mode. Different percentages of alpha-helices, beta-sheets, and TM-helices were found in the CmOMTs (Fig. S2). Nine CmOMTs (CmOMT7, CmOMT8, CmOMT10, CmOMT12–CmOMT15, CmOMT19, and CmOMT21) had the highest proportions of alpha helices (Fig. S2, Table S4). CmOMT2, CmOMT5 and CmOMT6 contained the highest number of beta-sheets, only CmOMT8, CmOMT9, CmOMT11, CmOMT14, CmOMT16, CmOMT18, and CmOMT21 contained TM-helixes (Fig. S2, Table S4). In addition, disordered regions (DRs) were found in all CmOMTs, and these regions warrant further exploration (Fig. S2, Table S4).

Synteny and GO enrichment analyses of CmOMT genes

Synteny analysis of CmOMTs was performed using the MCScan tool, and the fragment tandem duplication events of the CmOMT gene family members were investigated. As shown in Fig. S3A, only one CmOMT was involved in fragment tandem duplication events among the CmOMT gene family (Table S5). These results suggest that few CmOMTs arose from fragment tandem duplication events; thus, fragment tandem duplication events have a limited evolutionary role in CmOMTs. To further explore the phylogenetic mechanisms of CmOMTs, we constructed syntenic maps of melon based on three species, including a dicot (Arabidopsis thaliana) and two monocots (Citrullus lanatus and Cucumis sativus) (Fig. S3). Multispecies syntenic analysis showed that 2, 6, and 5 CmOMT genes had syntenic relationships with Arabidopsis, watermelon, and cucumber, respectively (Fig. S3). Collinear CmOMT9 and CmOMT14 gene pairs were found between melon and Arabidopsis (Fig. S3B, Table S5). Six collinear CmOMT gene pairs (CmOMT5/9/13/14/20/21) were identified between melons and watermelons (Fig. S3C, Table S5), five collinear CmOMT gene pairs (CmOMT5/9/14/18/21) were associated with cucumber (Fig. S3C, Table S5). Collinear CmOMT9 and CmOMT14 gene pairs were identified between melon and Arabidopsis/watermelon/cucumber, indicating that homologous gene pairs may have been generated before the divergence of monocotyledonous and dicotyledonous plants. Collinear CmOMT5 and CmOMT18 gene pairs were found between melon and other Cucurbitaceae species, suggesting that these OMT genes may have played key roles in the evolution of Cucurbitaceae OMTs.

To obtain a comprehensive understanding of the gene functions, GO term analyses were performed for all CmOMTs. The results of the GO enrichment analysis (Q-value ≤ 0.05) showed the genes corresponded to two biological process classes, namely, ‘methylation’ (Q-value = 6.92e−32) and ‘metabolic process’ (Q-value = 0.0055), and 11 molecular function classes (Table S6), which were closely associated with methyltransferase, such as ‘O-methyltransferase activity’ (Q-value = 2.6e−58), ‘methyltransferase activity’ (Q-value = 1.25e−33), and ‘caffeate O-methyltransferase activity’ (Q-value = 0.0016) (Table S6). These findings were consistent with previous reports showing that the OMT gene family is associated with methylation processes in plants (Liu et al., 2015).

Cis-acting element analysis of CmOMT genes in promoter regions

Cis-acting elements are transcription factor-specific binding sites that play important roles in regulating plant growth, differentiation, and development. We extracted a 2,000 bp sequence in the upstream region of each CmOMT and used the PlantCARE tool to identify the cis-acting elements in CmOMT promoters. Based on functional labeling, 10 cis-acting elements were obtained, among which 50% (5/10) was associated with hormonal responses and 50% was associated with stress and developmental responses (Fig. 3 and Table S7). Thirteen, ten, and seven CmOMT promoters contained ABA-responsive elements (ABREs), auxin-responsive elements (AuxRR-core and TGA-element), and GA-responsive elements (P-box, TATC-box, and GARE-motif), respectively (Fig. 3 and Table S7). MeJA-responsive elements (CGTCA-motif and TGACG-motif) and SA-responsive elements (TCA-element) were present in nine and eight CmOMTs promoters, respectively (Fig. 3 and Table S7). Wound-responsive elements (WUN-motifs) were present only in CmOMT2 and CmOMT8 promoters (Fig. 3 and Table S7). A drought-responsive element (MBS), defense and stress responsiveness element (TC-rich repeats), and low-temperature responsiveness (LTR) element were observed in 7, 8, and 4 CmOMT promoters (Fig. 3 and Table S7). The light-responsiveness element was commonly found in the promoters of all CmOMTs.

Expression analysis of CmOMT genes in different tissues and developmental stages

Our results revealed that a large number of development- and stress-related cis-acting elements were widely distributed in CmOMTs. We further explored the expression profiles of CmOMTs in different tissues and developmental stages. Gene expression profiling data were used to analyze the expression of CmOMTs in different tissues and at different developmental stages, and a gene expression heatmap was presented by scaling the expression of each gene across samples (Fig. 4), including the callus, dry seeds, root, stem (stem downside and upside), shoot apex, leaves (young, 6th, 9th and 12th leaves), tendril, flower (anther in male, petal in female, and stigma in female flowers), ovary (DAF 0, 2 and 4), fruit flesh (DAF 8, 15, 22, 29, 36, 43 and 50), and fruit epicarp (DAF 8, 15, 22, 29, 36, 43 and 50) (Table S8). Three CmOMTs (CmOMT16/18/21) were highly expressed in the callus, and only CmOMT14 was expressed in dry seeds. Nine CmOMTs were highly expressed in the roots, e.g., CmOMT5/18/21. Seven CmOMTs were highly expressed in the stems, with CmOMT1/10/18 showing high expression in both the downside and upside of stems. CmOMT9/10/14/18 was highly expressed in leaves (young, 6th, 9th, and 12th leaves), tendrils, and flowers (anther in male, petal in female and stigma in female flower). Notably, CmOMT9 was expressed at higher levels in the petals of female flowers than the other CmOMTs (Fig. 4 and Table S8). CmOMTs that play a key role in fruit pigment methylation can affect plant coloration and improve plant self-protection against environmental stress (Roldan et al., 2014; Du et al., 2015). In the ovary (DAF 0, 2, and 4), CmOMT14/18 was highly expressed. CmOMT1/18 showed an increasing and then decreasing trend during fruit flesh ripening (DAF 8, 15, 22, 29, 36, 43, and 50), consistent with the fruit epicarp. CmOMT14/15 showed an overall increasing trend in fruit flesh (DAF 8, 15, 22, 29, 36, 43, and 50), CmOMT14 displayed an increasing and then decreasing trend during fruit flesh ripening (Fig. 4 and Table S8). Overall, CmOMT1/14/15/18 were highly expressed throughout the growth and developmental stages of melon and in most tissues.

Expression of CmOMT genes in response to abiotic and biotic stresses

To investigate whether CmOMTs are involved in abiotic and biotic stress responses, we analyzed the expression patterns of 21 CmOMTs in varieties with different stress tolerances based on publicly available transcriptomic and qRT-PCR data. The stressors included salt (300 mM NaCl), HTH, cold (6 °C), Fusarium oxysporum f. sp. melonis race 1.2 (FOM1.2), and powdery mildew (Wang et al., 2016; Sebastiani et al., 2017; Zhu et al., 2018; Diao et al., 2020; Weng et al., 2022). The results showed that the expression levels of CmOMT6/9/12/14/18/21 were more variable after salt stress, which could be expanded on in a future study (Fig. 5A). The qRT-PCR results showed that CmOMT6 and CmOMT9 were initially downregulated and then upregulated within 7 d of salt stress, CmOMT12 was significantly downregulated, and CmOMT14/18/21 showed a peak in expression at 3 d after salt stress (Fig. 5D). In summary, CmOMT6/9/14/18 may have a positive regulatory role in response to salt stress, CmOMT14/18 may play a pivotal role in the early stages of salt stress, and CmOMT12 may have a negative regulatory role in response to salt stress. Based on transcriptome data, the results revealed that CmOMT14/18 was highly expressed after HTH stress (Fig. 5B). RT-qPCR analysis indicated that the expression of CmOMT14/18 was induced by HTH stress and that CmOMT14/18 were consistently upregulated after HTH stress (Fig. 5E). Of the 21 CmOMTs that responded to cold stress, four CmOMTs (CmOMT2/12/14/18) were highly expressed under both the control and cold stress conditions, suggesting that these genes may play a potential role in the mitigation of cold stress (Fig. 5D). RT-qPCR analysis revealed that the expression pattern of CmOMT2 was significantly decreased after cold stress and that the expression pattern of CmOMT14/18 was significantly upregulated after cold stress, indicating that these genes may have a positive regulatory role after cold stress (Fig. 5F). Overall, the transcript levels of CmOMT14/18 were significantly affected by multiple abiotic stressors, suggesting a positive regulatory role.

To explore the potential function of the 21 CmOMTs in response to biotic stress, published RNA-seq data were used to analyze the expression pattern of the CmOMT gene family under FOM1.2 and powdery mildew infection (Sebastiani et al., 2017; Zhu et al., 2018). Fifteen CmOMTs (CmOMT2/3/4/5/7/8/10/11/12/13/14/15/16/18/21) were significantly changed in ‘CHT’ (susceptible cultivar) after FOM1.2 infection, among which seven CmOMTs (CmOMT3/4/5/7/8/16/21) were significantly upregulated at 24 hpi (hours post FOM1.2 inoculation) and five (CmOMT3/5/14/16/18) were significantly upregulated at 48 hpi (Fig. 6A). Similarly, 15 CmOMTs (CmOMT2/3/4/5/6/7/8/9/10/12/13/15/16/18/20) were significantly up/downregulated in the resistant cultivar ‘NAD’ after FOM1.2 inoculation, only one (CmOMT18) was upregulated at 24 hpi, and eight (CmOMT3/4/5/7/8/13/16/20) were upregulated at 48 hpi (Fig. 6A). Interestingly, CmOMT18 reacted to FOM1.2 at an earlier stage in ‘NAD’ than ‘CHT,’ which showed changes at 48 hpi (Fig. 6A). Therefore, CmOMT18 may have a greater potential for FOM1.2 resistance. In addition, to verify whether the CmOMT gene family was involved in the response to FOM1.2, a correlation analysis was performed between resistance genes (PR-1 and LRR genes) and CmOMT s. PR-1 was positively correlated with CmOMT18, and most LRR genes showed a favorable positive correlation with CmOMT18 (Fig. S4). Taken together, CmOMT18 may play a positive regulatory role in the response to FOM1.2. To examine the reaction of the CmOMT gene family to powdery mildew, six CmOMTs (CmOMT3/4/5/7/11/18) were upregulated in MR-1 (resistant cultivar) after powdery mildew inoculation, with CmOMT3/7 upregulated at 1 and 3 dpi (days post-powdery mildew inoculation) and CmOMT5/18 upregulated at 3 and 7 dpi. Most CmOMTs were up/downregulated in the susceptible cultivar ’Top mark’ after powdery mildew inoculation: CmOMT2/9/10/12/15/16 were upregulated at 1 dpi, CmOMT3/9 was upregulated at 3 dpi, and CmOMT3 was upregulated at 7 dpi (Fig. 6B). As CmOMT18 was specifically expressed in leaves, a correlation analysis of the CmOMT gene family with R genes revealed that CmOMT18 may be involved in the regulation of melon resistance to powdery mildew (Fig. S4).