The GRAS gene family in watermelons: identification, characterization and expression analysis of different tissues and root-knot nematode infestations

Identification, characterization, and chromosomal location of watermelon ClGRAS genes

The newest versions of the genome annotations of watermelon were downloaded from the Cucurbit Genomics Database (http://www.cucurbitgenomics.org/). The latest HMM (Hidden Markov Model) model downloaded from the Pfam database (http://pfam.xfam.org/) was employed to search members of the GRAS transcription factor with a cut-off E-value of 1e⁻¹⁰ with GRAS domain (PFD03514) (El-Gebali et al., 2019). The GRAS protein sequences of Arabidopsis were obtained from the Arabidopsis Information Resource (https://www.arabidopsis.org/index.jsp) and rice GRAS sequences were obtained from Rice Gene Annotation Project (http://rice.plantbiology.msu.edu/index.shtml), and both of them were used as queries to explore the Cucurbit Genomics database by the default parameters. Combining with the utilization of the Pfam database and the Conserved Domain Database (CDD, https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi), only those sequences having a full-length GRAS domain were chosen as ClGRAS proteins and applied to the subsequent analyses.

The prediction of the physical and chemical features of ClGRAS proteins, such as molecular weight (MW) and the theoretical isoelectric point (PI) values were calculated by ExPASy (https://web.expasy.org/protparam/). To verify the subcellular localization of the identified ClGRAS proteins, LocTree3 Prediction system (https://rostlab.org/services/loctree3/) and WoLF PSORT Prediction (https://wolfpsort.hgc.jp/) were applied to predict the protein sequences. To identify the chromosomal locations of all ClGRAS genes in watermelon, the information of locus coordinates was obtained from the genomic sequences. MapChart software was used for the drafting of ClGRAS genes’ chromosomal positions and relative distances on the basis of their ascending order of physical position (bp) (Voorrips, 2002).

Multiple sequence alignment and phylogenetic analysis of GRAS proteins

Multiple sequence alignments of GRAS proteins in watermelon, rice, and Arabidopsis were performed using ClustalW (http://www.ebi.ac.uk/Tools/msa/clustalw2/). Arabidopsis and rice are generally used model plant species for exploring genetic relationships. MEGA X was further applied to establish an unrooted Maximum Likelihood (ML) phylogenetic tree with a pairwise deletion option and Poisson correction model (Kumar et al., 2018). Bootstrap analysis with 1,000 replicates was employed to examine the statistical reliability and choose the best results. The phylogenetic tree was visualized by the FigTree v1.4.3 program. Thereafter, members of watermelon ClGRAS gene families were further divided into diverse subfamilies based on well-established classification in other species (Huang et al., 2015; Liu & Widmer, 2014; Zhang et al., 2018).

Analysis of conserved motifs and gene structures

To determine the conserved motifs of each GRAS gene in watermelon, deduced GRAS protein sequences were submitted to MEME (version 5.1.1, http://meme-suite.org/tools/meme), with the default parameters except the number of motifs was chosen 6 (Bailey et al., 2009). To elaborate exon-intron organization for each ClGRAS gene, we downloaded the coding sequences (CDSs) and corresponding genomic sequences of ClGRAS genes in watermelon from the Cucurbit Genomics Database. Hereafter, both of them were subjected to the Gene Structure Display Server 2.0 (GSDS 2.0, http://gsds.cbi.pku.edu.cn), which was conducted and displayed by comparing CDSs and their corresponding genomic sequences.

Gene ontology annotation

With the intention for functional annotation identification of ClGRAS proteins, 37 watermelons ClGRAS proteins sequences were uploaded to Goplot (an R package) and were utilized for Gene Ontology (GO) analysis, and subsequently for mapping and annotation (Walter, Sanchez-Cabo & Ricote, 2015). The GO analysis of 37 watermelon ClGRAS proteins contains three parts such as biological process, molecular function, and cellular component.

Analysis of the protein–protein interaction network

Protein-protein interaction (PPI) data were acquired from the online database STRING (https://string-db.org/cgi/info.pl), an open-source software interface for predicting and visualizing sophisticated networks. The data for PPI was gathered from text mining of literature from peer-reviewed journals, containing the physical interactions and enzymatic reactions found in signal transduction pathways, computational predictions including those based on genomic context analysis as well as derived from analyses of co-expressed genes (Raman, 2010). The PPI data was preprocessed, containing removing redundancy and self-loops. Targets with a high confidence score >0.7 were selected to generate the PPI networks. PPI networks were visualized in Cytoscape software with the nodes representing proteins/genes and the edges indicating interactions between any two proteins/genes.

Expression analysis of ClGRAS genes during watermelon fruit development and root-knot nematode infection

The transcriptome data of watermelon fruit flesh and rind tissue sampled at 10, 18, 26, and 34 days after pollination (DAP) under BioProject ID SRP012849, were downloaded from the National Center for Biotechnology Information Sequence Read Archive (Bethesda, MD, USA) database. The reference transcriptome was chosen all the annotated gene sequences of the reference genome from watermelon cultivar 97103. In addition, for assessing the expression levels of watermelon ClGRAS genes during RKN infection, the transcriptome data of RKN infection were downloaded from the genome sequence Archive in the BIG Data Center GSA database (https://bigd.big.ac.cn/gsa/browse) under the accession numbers of CRA001311 and CRA001312. It was composed of four different treatments, CK (white light and water solution), RL (red light treatment and water solution), RKN (white light and M. incognita infection), and RL+RKN (RR, red light treatment, and root-knot nematode M. incognita infestation). Gene expression levels were estimated and normalized with TopHat/Cufflinks pipeline with FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values extracted from the above-mentioned transcriptome data (He et al., 2015). The log2-transformed FPKM values were employed to draw a heatmap to depict the expression of each ClGRAS gene by using the RStudio (version 1.1.463; RStudio, Boston, MA, USA).

To further comprehensively analyze the expression profiles of ClGRAS genes in watermelon, we performed qRT-PCR experiments of four randomly selected genes (including ClGRAS2, 18, 28, and 33). We extracted total RNA from roots and leaves under four treatments (CK, RKN, RL, and RR, respectively) via using the TRIzol Reagent (Takara, Beijing, China) following the manufacturers’ recommendations. Then, Reverse transcription was performed by PrimeScript RT reagent Kit (TaKaRa, Beijing, China). CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) was used to perform the qRT-PCR analysis. The watermelon 18S ribosomal RNA (18S rRNA) was chosen as the endogenous control (Dong & Shang, 2013). Specific primers of ClGRAS genes used for qRT-PCR were designed by Primer Premier 6.0 and listed in Table S1. Specifically, qPCR was performed in a 15 μL reaction mixtures containing 7.5 μL 2× SYBR® Premix Ex Taq™ II (TaKaRa, Beijing, China), 1 μL cDNA template, 0.3 μL each gene-specific primer, and 5.9 μL ddH2O. The reaction conditions were 30 s at 94 °C, 45 cycles of 20 s at 94 °C, 20 s at 55 °C, and 30 s at 72 °C. The melting curves were analyzed from 60 °C to 95 °C to observe the specificity of the PCR products. Three biological replicates were conducted for each treatment, and each composed of three technical replicates. Expression levels were determined as the mean across the three replicates and the relative expression between samples was calculated using the comparative 2^−ΔΔCT method (Livak & Schmittgen, 2001). The relative expression levels of CK were determined as 1 for comparison.

Genome‑wide identification and characterization of GRAS family members in watermelon

Candidate ClGRAS genes were identified from the watermelon genome using the Base Local Alignment Search Tool (BLAST) with the query sequences of Arabidopsis and rice GRAS genes (Altschul et al., 1997). Subsequently, the retrieved sequences were submitted to the CDD and Pfam databases to confirm the presence of conserved domains. A total of 37 candidate ClGRAS genes were identified in the watermelon genomes and named according to their chromosomal locations (Table 1). The length of 37 ClGRAS proteins ranged from 226 (ClGRAS32) to 859 (ClGRAS15) amino acid residues.

The distribution of 37 ClGRAS genes on 11 chromosomes was visualized by the MapChart program. It can be seen that all of the 37 ClGRAS genes are unevenly distributed. In addition, most of ClGRAS genes appeared to congregate on the proximate or the distal ends of the chromosomes (Fig. 1). There were 7, 5, and 5 ClGRAS genes mapped onto chromosomes 1,9 and10 respectively, whereas chromosomes 4, 8, and 11 harbored just one ClGRAS gene.

The ExPASy online tool was used to predict the basic physicochemical properties of 37 ClGRAS proteins. As the results showed in Table 1, the molecular weights of these ClGRAS proteins were ranged from 31.08 (ClGRAS32) to 92.89 (ClGRAS15) kDa. Moreover, the majority of the ClGRAS proteins exhibited acidic isoelectric points (less than 6.94), with the lowest being 4.50 (ClGRAS6), only two proteins had alkaline isoelectric points of more than 7.47 (ClGRAS31), of which ClGRAS32 was the highest at 10.16.

Based on the predicted subcellular localization results, ClGRAS proteins were located in the nucleus, chloroplast, and cytoplasm, more than half of it distributed in the nucleus, illustrating that these proteins as transcription factors play a transcriptional regulatory role directly in the nucleus. Detailed information of ClGRAS genes in the watermelon is presented in Table 1.

Classification and phylogenetic analysis of the GRAS genes

To evaluate the phylogenetic relationships among the ClGRAS gene members, all the genes from watermelon, rice, and Arabidopsis were aligned by the Maximum Likelihood method to construct an un-rooted phylogenetic tree. ClGRAS genes in the three species were classified into ten groups as shown in the phylogenetic tree in Fig. 2. Eight watermelon ClGRAS proteins were belonging to the HAM subfamily and constituted the largest subgroup of the phylogenetic analysis, while the DTL subfamily comprised only one member (ClGRAS33). Concerning the HAM subfamily, ClGRAS genes were clustered and homologous with rice and Arabidopsis GRAS genes, implying that they could be derived from gene duplications of the same gene locus during the evolution of the watermelon genome. These results indicated that these paralogous genes may play an identical functional role in the growth of watermelon. Apart from this, the un-rooted tree containing 6, 5, 3, 2, 4, 3, 2, and 3 ClGRAS members in the SHR, PAT1, LISCL, SCL3, DELLA, SCR, Os4, and SCL4/7 subcategories, respectively (Fig. 2).

Motif composition and Exon-intron structures of ClGRAS gene families

A total of six conserved motifs were identified in the ClGRAS proteins and their detailed consensus sequence information is listed in Table 2. According to the phylogenetic tree and conserved motifs (Fig. 3A), we could know that the same group of ClGRAS genes shares similar conserved motifs organization concerning either motif number or gene length, which revealed that there could be similar genetic functions. For instance, ClGRAS3, ClGRAS18, and ClGRAS26, ClGRAS27, and ClGRAS15 have 5 motifs within the same subgroup. There were 4 motifs (motif 2, 3, 4, 6) in the C-terminal region of most ClGRAS proteins, while motif 1 and 5 were found in the N-terminus. Almost all ClGRAS proteins harbored motif 2, 3, 4, and 6 in the C-terminal region, while the motifs across the different proteins were variable in N-terminal. These results are in line with previous research that the C-terminus might be more conserved than the N-terminus via evolution (Fuxreiter, Simon & Bondos, 2011). Moreover, it has previously been demonstrated that variety in C-terminal domains of GRAS proteins can increase the functional versatility and the complexity of biological networks of the GRAS proteins (Sun, Jones & Rikkerink, 2012; Sun et al., 2011). Furthermore, the amino acid frequency at six conserved motifs of ClGRAS proteins was investigated and then downloaded the LOGO of these protein motifs. It indicated that the amino acid frequency of the conserved motifs is not consistent in different ClGRAS proteins (Fig. S1).

In order to verify the reliability of the phylogeny analysis and to elucidate the protein functions, as well as structural diversity of ClGRAS genes, gene structures, and conserved motifs in the 37 ClGRAS proteins, were identified using the GSDS 2.0 website and MEME program. As shown in Fig. 3B, the number of introns ranged from zero to 2. A considerable proportion of ClGRAS genes do not contain introns, and an analogous phenomenon of deletion of introns was observed in the genomes of other species. Eight genes such as ClGRAS9, ClGRAS13, ClGRAS15, ClGRAS19, ClGRAS20, ClGRAS31, and ClGRAS37 possess one intron each, whereas just one gene (ClGRAS11) consists of 2 introns. It should be noted that some genes such as ClGRAS13 and ClGRAS31shared similar intron numbers but with different intron lengths. On the whole, ClGRAS members belonging to the identical subfamily on watermelon had similar exon-intron structures.

Gene ontology annotation

In order to gain further insight into which biological processes these ClGRAS genes involved, a GO annotation analysis of the ClGRAS protein was performed and the GO number was present in Table S2. According to the results obtained, ClGRAS proteins were involved in many biological processes (Fig. 4). The many ClGRAS proteins are responsible for cellular processes (GO: 0009987), regulation of biological processes (GO: 0050789), and biological regulations (GO: 0065007). Also over 10 proteins were involved in the metabolic processes (GO: 0008152). It is worth noting that the number of ClGRAS proteins localized in the organelle (GO: 0043226) is the same as that of the cell (GO: 0044464) and cell part (GO: 0005623), which showed that ClGRAS proteins perform a functional role in these locations. Additionally, a moderate number of ClGRAS proteins might be associated with transformation regulatory activity (GO: 0140110). Nonetheless, just a few ClGRAS proteins respond to stimulus (GO: 0050896) and play role in signaling (GO: 0023052) and reproductive processes (GO: 0022414). Thus, according to the above-mentioned analysis of the biological processes of ClGRAS genes, the multi-functional role of watermelon ClGRAS proteins could be involved in various biological processes, such as cellular/metabolic processes, and had functions in processes of signal transduction as well as the stress response.

Interaction network of ClGRAS proteins in watermelon

To investigate the potential molecular mechanisms of watermelon ClGRAS proteins, the data only experimentally determined in the STRING database were selected to construct the PPI network. By using Cytoscape software, we found that a total of 25 interaction protein pairs were predicted and nine ClGRAS proteins were included in this network (Fig. 5). Apart from these ClGRAS proteins, eight proteins including PHYA and PHYB (Phytochrome A and B), JKD (C2H2-like zinc finger protein), GID1B and GID1C (alpha/beta-Hydrolases superfamily protein), BZR1 (Encodes a positive regulator of the brassinosteroid signalling pathway), HD1 (Encodes a histone deacetylase) and DAG1 (Dof-type zinc finger DNA-binding family protein) were also presented in this network. Among them, ClGRAS16, JKD, and BZR1 were paired with 4 proteins, separately; while ClGRAS11, ClGRAS12, GID1B, and GID1C were associated with three proteins, and the remaining proteins interacted with other one or two proteins. Apparently, two ClGRAS proteins exhibited the most protein-protein interactions, ClGRAS1 and ClGRAS29, which revealed that these two proteins played a key role in the whole protein interactions network. Moreover, ClGRAS16, JKD, and BZR1 featured prominently in the protein-protein network. It was illustrated that the ClGRAS protein correlated with them may possess the function of a nuclear-localized putative transcription factor with three zinc finger domains (JKD) and regulating the brassinosteroid signalling pathway (BZR1). In summary, the results obtained in the present study showed various types of protein interactions of the ClGRAS proteins. Although the ClGRASs’ interaction network needs to be further verified by experiment, our results will further expedite the identification of more essential proteins and biological modules which interacted with ClGRAS proteins and provide vital theoretical evidence for the molecular mechanisms of ClGRASs. The detailed information of the proteins in this PPI network is presented in Table S3.

Expression analysis of ClGRAS genes in fruit flesh and rind during watermelon fruit development

To further evaluate the potential functions of ClGRAS genes in the growth and development of watermelon, transcriptome data derived from the fleshes and rinds of watermelon fruit at four developmental stages were applied to determine the expression levels of ClGRAS gene families. From the heat map drawn by the Rstudio (Fig. 6), a total of 26 ClGRAS genes were not expressed during the development of watermelon fruit, which might be lost their functions during evolution. Furthermore, the remaining 11 ClGRAS genes certainly participated in the regulation of fruit development and maturity. Overall, the expression patterns of these 11 genes were altered, either up-regulated or down-regulated. Of them, ClGRAS7, ClGRAS8, ClGRAS16, and ClGRAS26 had relatively high expression in the fruit flesh and rind at 10 DAP and 18 DAP, respectively, while relatively low expression levels in fruit flesh at 26 and 34 DAP. Apart from that, there existed 5 ClGRAS genes which displayed a relatively low-level expression in fruit flesh at 26 and 34 DAP. Additionally, ClGRAS1 and ClGRAS17, ClGRAS24, and ClGRAS28 were significantly expressed in fruit flesh at 26 DAP and in the rind at 10 DAP, respectively, which indicated that these genes may function in rind and fruit flesh development of watermelon. According to Gene Ontology, these differential expression genes are annotated with the 5 to 17 GO terms, such as transcription, DNA-templated (GO: 0006351), regulation of transcription, DNA-templated (GO: 0006355), and biological_process (GO: 0008150). Surprisingly, three genes, ClGRAS2, ClGRAS19, and ClGRAS25 exhibited a similar expression in the rind at 18, 26, and 34 DAP, suggesting that these three genes may have a synergistic effect in the rind development. Also, several genes have no expression at some time points, such as ClGRAS1, which was not expressed at rind tissue at 18 DAP. Collectively, overlapping but diverse expression levels of the ClGRAS genes showed that these genes might perform a variety of crucial roles in the growth and developmental processes of watermelon.

Expression analysis of ClGRAS genes under root-knot nematode infection

To gain insight into the potential functions of ClGRAS genes in response to biotic stresses, expression profiles of each ClGRAS gene in the roots and leaves of watermelon under root-knot nematode infection were investigated. As shown in Fig. 7A, under the treatments of CK, RKN, RL, and RR, expression levels of the ClGRAS genes were varied in roots. The expression of 11 ClGRAS genes (2 up-regulated and 9 down-regulated) was significantly changed by RKN treatment compared with CK. In addition, a total of 10 and 3 ClGRAS genes were respectively up-regulated and down-regulated by RL treatment compared with RKN. Moreover, compared with RL treatment, 3 and 8 ClGRAS genes were found to be (up-regulated and down-regulated, respectively) altered expression under RR treatment. Among the expression profiles of these ClGRAS genes, it was obviously found that 10 genes had high expression under RL treatment, whereas reduced in other treatment. These results illustrated that the genes played a vital role in response to red light treatment. The analogous phenomenon has existed in other genes, for instance, ClGRAS18 with RR and ClGRAS8 on RKN treatment.

We also determined the expression levels of the ClGRAS genes in leaves, and the results of changes in expression profiles are shown in Fig. 7B. A total of 8 and 3 ClGRAS genes showed apparent upregulation and downregulation under the treatment of RKN relative to CK, respectively. Additionally, there were 12 and 10 ClGRAS genes with remarkably incrementing and decreasing expression by RL treatment compared with RKN. In contrast, under the treatment of RR, only five genes had a relatively high expression, the rest of the genes were either only expressed at very low levels or not expressed. Similar to the root, there were a string of genes that were highly expressed under RL treatment, while were not expressed or only expressed at very low levels under other treatments, which revealed that these genes are putatively involved in the response to RL treatment. Similar results were also found in RKN treatment. Furthermore, based on Gene Ontology, these high expression genes are annotated with the 8 to 19 GO terms, for example, transcription, DNA-templated (GO: 0006351), transcription factor activity, sequence-specific DNA binding (GO: 0003700), cellular_component (GO: 0005575), molecular_function (GO: 0003674), sequence-specific DNA binding (GO: 0043565), regulation of transcription, DNA-templated (GO: 0006355), and biological_process (GO: 0008150).

To test the accuracy of gene expression determined from transcriptome data, we analyzed the expressions of 4 randomly selected ClGRAS genes under four different treatments by qRT-PCR (CK, RKN, RL, and RR). As showing in Fig. 8, expression levels of the 4 ClGRAS genes differed among four treatments, indicating their involvement in response to these treatments. In roots, relative expression levels of ClGRAS2 and ClGRAS18 were differentially down-regulated under another three treatments compared with CK. ClGRAS28 exhibited closely resemble expression levels by CK, RKN, and RL treatments, while showed remarkably high expression levels by RR treatment. Noticeably, the relative expression level of ClGRAS33 was extraordinarily altered by RKN, RL, and RR treatments, especially by RR treatment. Moreover, expression profiles were changed in leaves, ClGRAS2 was significantly upregulated after RKN, RL, and RR treatments, whereas down-regulated by RR treatment when compared with RKN and RL (Fig. 9). ClGRAS18 held similar expression levels by RKN and RL treatments and evidently restrained by RR treatment. It is noteworthy that ClGRAS28 and ClGRAS33 have analogous expression patterns and higher expression levels by RKN and RR treatments, however, lower in RL treatment. On the whole, the qRT-PCR results were in accord with that of previous transcriptome results.