Identification and expression profile analysis of the SnRK2 gene family in cucumber

Plant materials

Cucumber seeds L306 (Cucumis sativus L, “L306” cultivar) were used in this experiment. The seeds were disinfected with 5% sodium hypochlorite and then germinated until the two true leaves were completely expanded. The seedlings were cultured in Yamazaki cucumber nutrient solution at 25 °C/18 °C (day/night), 250 µmol.m⁻² s⁻¹ light and 70% relative humidity in an artificial climate chamber. Samples were collected at three true leaf stages. This experiment set up four treatments: 50 µM ABA, 100 µM ABA, 200 mM NaCl and 10% PEG. In each treatment, a total of 20 seedling samples of similar size were treated. Five sampling time points were set for different treatments, and the leaves of four seedlings were collected at each time point (0, 3, 6, 12, and 24 h). RNA was extracted after cutting and mixing, and three CT values were measured as repetition.

Identification of CsSnRK2 genes

To screen candidate SnRK2 genes in cucumber, BLASTp was used to compare the protein sequences of 10 SnRK2s of Arabidopsis with that of cucumbers using an online database (http://cucurbitgenomics.org/). Set the BLASTp e-value to 1e⁻¹⁰ (Altschul, 1997). The redundancy of the results was removed manually. The hidden Markov model (HMM) of the kinase domain (PF00069.25) was obtained from the Pfam database (http://pfam.xfam.org/), and the cucumber genome database was further searched by HMMER 3.0 software. TBtools was also used to screen candidate genes; Pfam (http://pfam.xfam.org/search#tabview=tab1) and NCBI-CDD (https://www.ncbi.nlm.nih.gov/cdd/) databases were used for validation and candidate genes without SnRK2 domains (registration number: PF00069.25) were manually deleted (Zhou et al., 2018).

Chromosome localization and basic information analysis

For the distribution on chromosomes, the gff3 file of ”9930” V3 version of China long was downloaded. SnRK2 genome location information was obtained by screening. It was named CsSnRK2 s according to its distribution on chromosomes. The chromosome map of the gene was drawn by MG2C (http://mg2c.iask.in/mg2c_v2.0/). WoLF PSORT (https://wolfpsort.hgc.jp/) was used to predict protein subcellular localization (Xiong et al., 2015). The number of amino acids as well as the molecular weights, isoelectric points, and other physical and chemical information of CsSnRK2 gene were analyzed on ExPaSy website (https://web.expasy.org/protparam/) (Liu et al., 2017). The secondary structure of the cucumber SnRK2 protein was predicted by using the PRABI website (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=npsa_sopma. html).

Collinearity analysis of CsSnRK2 genes

Collinearity analysis using MCScanx provided by TBtools software. Cucumber genome was compared with Arabidopsis genome, maize genome and rice genome by BLASTp, and the search threshold was set as e-value <1e⁻⁵. Other default parameters have not been modified. The results of genome-wide BLASTp were used to calculate all possible collinear pairs between chromosomes. Using TBtools to map the collinear pairs of SnRK2 genes in cucumber, Arabidopsis, maize and rice (Chen et al., 2020). The protein sequences of the gene pairs were extracted for multiple sequence alignment, and then the alignment results and the CDS sequences of the gene pairs were used to calculate the non-synonymous/synonymous (d_N/d_S) values of duplicate gene pairs using PAL2NAL (http://www.bork.embl.de/pal2nal/index.cgi?) (Wang et al., 2021).

Phylogenetic, structure and conserved motif analysis

The phylogenetic trees of SnRK2 gene family of cucumber, A.thaliana, rice, maize, sugar beet, tomato, grapevine and pak-choi were constructed by using MEGA 11 software with ClustalW methods to align multiple sequence (Kumar, Stecher & Tamura, 2016). The bootstrap values of 1000 replicates were calculated at each node. The minimum neighbor method was selected to represent the phylogenetic tree. Then, the phylogenetic tree was visualized using the iTOL online tool (http://itol.embl.de/) (Ivica & Peer, 2016). The CDs and whole genome sequences of 12 cucumber SnRK2 genes were uploaded to GSDS (http://gsds.cbi.pku.edu.cn/) to visualize the positions of introns and exons and untranslated regions (Hu et al., 2014). We analyzed the conserved motifs of gene by MEME (http://meme-suite.org/tools/meme). The maximum motif number was set to 10, and the other parameters are the default values (Bailey et al., 2009).

Analysis of cis-acting elements, tissue-specific expression and protein interaction networks

TBtools were used to extract the genomic sequence 2,000 bp upstream of the start codon (ATG) of SnRK2 genes in cucumber. Then, the 2000 bp sequences of 12 CsSnRK2 were submitted to PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) for prediction. The RNA-Seq data of cucumber in different tissues and organs (tendril-base, tendril, root, leaf, stems, ovary-unfertilized, ovary-fertilized, ovary) were obtained from the cucumber genome data using the registration number PRJNA80169 (http://cucurbitgenomics.org/). The data were converted by log2 method, and TBtools software was used to draw the heat map of CsSnRK2 genes expression. Using 12 CsSnRK2 protein sequences as targets, the protein-protein interaction network was predicted by using STRING website (https://string-db.org/cgi).

RNA extraction and Real-Time PCR (qRT-PCR)

Total RNA was isolated from collected samples using a Plant RNA Extraction kit (Tiangen, China). cDNA was synthesized with fastking cDNA dispersion RT supermaxs Kit (Tiangen, China) and 2uL RNA was used as template. Design of primers using CDs sequences of genes. The SYBR Green kit (Tiangen, China) was used for fluorescence quantitative analysis. The volume of reaction system was 20 uL, which contained 2 uL cDNA solution, 10 uL 2*SuperReal PreMix Plus, 0.6 uL of 10 uM forward and reverse primers, 0.4 uL 50*ROX Reference Dye and 6.4 uL of distilled deionized water. qRT-PCR analysis using LightCycler^® 480 II real-time fluorescent quantitative PCR instrument. Each qRT-PCR reaction was carried out with three technical replicates. Cucumber actin (DQ115883) was used as internal control. The amplification primers and internal reference primers are shown in Table 1. Amplification conditions were as follows: 95 °C for 15 min, and 40 cycles of 95 °C for 10 s and 60 °C for 30 s. The relative expressions levels of the SnRK2 genes were calculated using the 2^−ΔΔCt method (Livak & Schmittgen, 2002). SPSS 20.0 software was used for one-way ANOVA, and Duncan method was used for significance test (p < 0.05). Using origin 9.0 software to complete the histogram of relevant expression, and the values of the last drawn graph are the average values of three repetitions.

Chromosome distribution and basic information of CsSnRK2 genes

To explore the function of the SnRK2 gene, we systematically analyzed the physicochemical properties of the cucumber SnRK2 protein. We identified 12 CsSnRK2 genes and named them sequentially from CsSnRK2.1 to CsSnRK2.12. The SnRK2 genes differed substantially according to the encoded protein size and biophysical properties (Table 2). The number of amino acids ranged from 298 to 673. The sequence of CsSnRK2.6 protein was the longest, containing 673 amino acids, while the sequence of CsSnRK 2.9 protein was the shortest, containing 298 amino acids. The molecular weight of CsSnRK2 ranged from 34,754.37 to 77,009.99 Da, and the theoretical isoelectric point (pI) ranged from 4.71 to 8.8. CsSnRK 2.5 and CsSnRK 2.7 were basic proteins, whereas the others were acidic proteins. The total average hydrophobic index of the 12 CsSnRK2 gene family members was negative, indicating that they were hydrophilic proteins.

Subcellular localization prediction (Table 3) results showed that all members of the cucumber SnRK2 gene family were localized in the cytoplasm, and most were also localized in the nucleus. The predicted secondary structure of the cucumber SnRK2 protein was mainly composed of an α-helix (36.05%–44.97%) and irregular curl (36.10%–45.43%), and the extended structure (9.40%–20.99%) and ß-turn (3.76%−7.37%) were relatively small.

Chromosome analysis of the cucumber SnRK2 gene family showed that 12 CsSnRK2 genes were unevenly distributed on five cucumber chromosomes (Fig. 1). There were two genes distributed on chromosome 2 (CsSnRK2.1 and CsSnRK2.2) and one gene each on chromosomes 3 (CsSnRK2.3), 6 (CsSnRK2.11), and 7 (CsSnRK2.12). The other genes were located on chromosome 4 (CsSnRK2.4-CsSnRK2.10). No members of this gene family were found on chromosomes 1 or 5.

Collinearity and selective pressure analysis of the CsSnRK2 genes

We explored the evolutionary process of CsSnRK2, AtSnRK2, ZmSnRK2, and OsSAPK2 genes. Collinearity analysis using MCScanx was performed using the TBtools software (Fig. 2). There were three pairs of collinearity genes between cucumber and A.thaliana, four pairs between cucumber and rice, and one pair between cucumber and maize. The col-linearity of the AtSnRK2, OsSAPK, ZmSnRK2, and CsSnRK2 genes could be broadly classified into two categories; one-to-many (CsSnRK2.11-AtSnRK2.7/8, CsSnRK2.11-OsSAPK1/2), or one-to-one (CsSnRK2.2-AtSnRK2.2, CsSnRK2.12-OsSAPK8, CsSnRK2.3-OsSAPK3, CsSnRK2.12- ZmSnRK2.8). In addition, these gene pairs belonged to the same group in the phylogenetic tree (Fig. 3), and their gene structure and conserved motifs were highly similar (Fig. 4).

Purifying selection is the process of removing unfavorable mutations, whereas positive selection is the accumulation of new favorable mutations and their transmission to the entire population (Qiao et al., 2015). The value of d_N/d_S plays a crucial role in gene selection and evolution; d_N/d_S >1 is positive selection, d_N/d_S = 1 is neutral selection, 0 < d_N/d_S <1 is purifying selection (Yadav et al., 2015). The results showed that all three pairs of genes had d_N/d_S values <1; indicating that the evolutionary process underwent purifying selection pressure and evolved conservatively, which helped maintain the basic functions of the genes (Table 4).

Phylogenetic tree analysis of the CsSnRK2 genes

To explore the functional characteristics and evolutionary relationships of the CsSnRK2 genes, a phylogenetic tree was constructed using multiple sequences alignments of SnRK2 protein sequences from the 12 cucumber, 10 Arabidopsis, 10 rice, 11 maize, six sugar beet, seven tomato, eight grapevine, and 13 pak-choi genes. The phylogenetic tree showed that the SnRK2 gene family was clustered into three branches, labeled as group I, group II, and group III (Fig. 3), which is consistent with the results of previous studies in Arabidopsis (Umezawa et al., 2009). The cucumber SnRK2 family members of cucumbers were unevenly distributed among the three groups. Group I had eight members, group II had one member, and group III had three members. The phylogenetic tree showed that CsSnRK2 was more closely grouped with SnRK2 from the dicots Arabidopsis, sugar beet, tomato, grapevine, and pak-choi than with SnRK2 from the monocots rice and maize. The cucumber SnRK2 and grapevine SnRK2 proteins were highly homologous, indicating that they have a close phylogenetic relationship and conservative evolution. Specifically, CsSnRK2.3 and VvSnRK2.6; CsSnRK2.11 and VvSnRK2.4/8; and CsSnRK2.2 and VvSnRK2.1 had the closest phylogenetic relationships, indicating that they may have the most similar functional characteristics.

Structure and conserved motif analysis of the CsSnRK2 genes

Analysis of the gene structure and conserved motifs of cucumber SnRK2 provided important information for understanding the evolution of this gene (Fig. 4). Therefore, a phylogenetic tree was constructed using the CsSnRK2 and AtSnRK2 genes. According to the phylogenetic tree, the cucumber and Arabidopsis SnRK2 gene families could be divided into three groups (Fig. 4A). In the third group, all cucumber genes differed from the previous (Fig. 3). The exon-intron distribution of cucumber SnRK2 genes was conserved, most of these had nine exons and eight introns; this was consistent with most Arabidopsis SnRK2 genes, except for AtSnRK2.6 (10 exons), AtSnRK2.3 (10 exons), and AtSnRK2.8 (6 exons) (Fig. 4B) (Yoon et al., 1997). In addition, the number and location of introns and exons in the third group of CsSnRK2 were different from those in Arabidopsis homologs, which may contribute to the functional diversification of CsSnRK2. Moreover, most introns of the cucumber SnRK2 gene were found to be much longer than those in the Arabidopsis SnRK2 genes.

We used TBtools to predict the conserved motifs in the CsSnRK2 and AtSnRK2 families (Fig. 3C). A total of 10 conserved motifs were identified, and their motif sequences are shown in Fig. 4D. Except for AtSnRK2.3, all genes had motifs 1-10, and the sequence and position of the motif in the structure were nearly the same, showing high conservatism in evolution. In the third group, there were significant differences in motif placement and location, and motif deletions were also observed. For example, CsSnRK2.7 and CsSnRK2.8 lacked motif 9, while CsSnRK2.9 lacked motif 4/5/6 but hads two of motif 1, which was similar to that of CsSnRK2.5. We hypothesized that these differences were associated with gene-specific functions.

Cis-acting element of the CsSnRK2 promoter region

The plant SnRK2 gene family plays an important role in abiotic stress responses, and the analysis of regulatory elements provides favorable information for the study of gene function. Several cis elements were identified in the promoter regions of the CsSnRK2 genes (Fig. 5). Among them, at large variety of hormone-related cis elements, such as P-box, TATC-box, GARE-motif, TCA-element, ABRE, TGACG-motif, CGTCA-motif, TGA-element, and AuxRR-core, indicated that they responded to different hormones such as MeJA, growth hormone, ABA, ethylene, GA and SA. The four cis elements were MBS, LTR, TC-rich repeats and ARE, which are involved in drought, low temperature, defense, stress and anaerobic induction. In addition, six cis elements associated with development were identified in the CsSnRK2 promoter, including P-box, TATC-box, Gare-Motif, ABRE, and TGA-element/AuxRR-core. The CsSnRK2.8 gene had the least active elements, containing only three elements, while the other genes had at least five elements. The genes, except CsSnRK2.7 and CsSnRK2.8, contained ABRE elements. This suggests that most CsSnRK2 genes are responsive to ABA. Only four members (CsSnRK2.1, CsSnRK2.2, CsSnRK2.4, and CsSnRK2.7) contained the LTR-element. In addition, two members (CsSnRK2.9, and CsSnRK2.11) contained TC-rich repeats, and five members (CsSnRK2.3, CsSnRK2.7, CsSnRK2.8, CsSnRK2.11, and CsSnRK2.12) contained MBS. We hypothesized that exogenous stress could induce CsSnRK2 gene expression through cis-responsive elements, thereby enhancing plant resistance to adversity.

Tissue-specific expression analysis of the CsSnRK2 genes

The heat map of RNA-Seq data expression of the CsSnRK2 genes in different tissues and organs using TBtools software (Fig. 6) revealed that CsSnRK2.1 and CsSnRK2.6 were not expressed in almost all tissues. CsSnRK2.2, CsSnRK2.4, CsSnRK2.8, CsSnRK2.10 and CsSnRK2.12 were strongly expressed in each tissue, and their function may be related to growth and development. CsSnRK2.5 was strongly expressed in the roots and moderately expressed in unfertilized ovaries. CsSnRK2.7 was strongly expressed in the roots and weakly expressed in other tissues. It was speculated that CsSnRK2.5 and CsSnRK2.7 may play an important role in root growth. CsSnRK2.11 was highly expressed in tendrils and the base of tendrils, and moderately expressed in other tissues. CsSnRK2.3 was moderately expressed in all the tissues. CsSnRK2.9 was moderately expressed in the roots and weakly expressed in the leaves.

qRT-PCR expression characteristics of the CsSnRK2 genes

The expression levels of the CsSnRK2 genes under salt, PEG, and ABA stress were determined using qRT-PCR. A total of 50 µmol/L ABA (T1), 100 µmol/L ABA (T2), 200 mmol/L NaCl (T3), and 10% PEG (T4) were sampled at 3, 6, 12, and 24 h, respectively, and 0 h was used as a control. The expression of CsSnRK2 was significantly different between the different time points and treatments (Fig. 7).

Under 50 µmol/L ABA treatment (Fig. 7A), four genes (CsSnRK2.4/5/11/12) were up-regulated, of which CsSnRK2.11, and CsSnRK2.12 were up-regulated the most, and their relative expression levels at 6 and 24 h were 4 and 3.5 times higher than that at 0 h, respectively. CsSnRK2.4 was highly expressed at 24 h, and its relative expression level at 24 h was 4.3 times higher than that at 0 h. In addition, five genes (CsSnRK2.2/6/8/9/10) were down-regulated, and CsSnRK2.1/3/7 were up-regulated first and then down-regulated. CsSnRK2.3/7 genes reached their maximum expression values at 3 h.

Upon 100 µmol/L ABA treatment (Fig. 7B), all the genes were up-regulated and then down-regulated, except for CsSnRK2.2 and CsSnRK2.5, which were downregulated. All three genes (CsSnRK2.4/11/12) were up-regulated, of which CsSnRK2.11 was up-regulated the most, and the relative expression level at 24 h was 11.9 times higher than that at 0 h. The expression of CsSnRK2.1 was the highest at 6 h, and the relative expression levels at 6 h was 6.8 times higher than that at 0 h. The maximum expression of CsSnRK2.2/6/7/9/10 was reached after 3 h; however, the expression of CsSnRK2.2 was down-regulated by both ABA treatments, suggesting that CsSnRK2.2 may play a negative regulatory role.

In response to 200 mmol/L NaCl treatment (Fig. 7C), three genes (CsSnRK2.4/11/12) were up-regulated, of which CsSnRK2.11 and CsSnRK2.12 were up-regulated the most; their relative expression levels at 12 and 24 h were 7.3 and 5 times higher than that at 0 h, respectively. In addition, five genes (CsSnRK2.2/6/8/9/10) were down regulated. The genes CsSnRK2.1 and CsSnRK2.7 were up-regulated first and then down-regulated.

For the 10% PEG treatment (Fig. 7D), six genes (CsSnRK2.1/5/7/10/11/12) were up-regulated, where CsSnRK2.11 and CsSnRK2.12 were up-regulated the most; their relative expression level at 12 and 24 h were 15.7 and 9.6 times higher, respectively, that of 0 h. However, two genes (CsSnRK2.2/4) were also down-regulated. The expression of CsSnRK2.5 reached a maximum at 12 h, and its relative expression level at 12 h was 6.6 times higher than that at 0 h.

This study provides an overview of CsSnRK2 expression under ABA, NaCl, and PEG stress conditions, where different CsSnRK2 members had different expression patterns under these stressors. This illustrates that the functional diversity of this gene is widely present in different plant species.

Analysis of the protein interaction network of the CsSnRK2 family genes

A protein–protein interaction network between CsSnRK2 and other cucumber proteins was constructed using the STRING database. According to the prediction results, cucumber SnRK2 protein can interact with five other proteins (XP_004146814.1, XP_004148120.1, XP_004167207.1, XP_004154794.1, and XP_004168589.1), the first four of which are PP2C proteins, which belong to the PP2C superfamily, and the fifth is a sucrose non-fermentable protein, presumably belonging to the other subfamilies outside SnRK2. The results also showed that eight CsSnRK2 proteins interacted with XP_004146814.1 (PP2C2), XP_004148120.1, and XP_004167207.1. The XP_004154794.1 protein interacted with five CsSnRK2 proteins. In addition, no interaction between CsSnRK2 proteins were found in the prediction results; however, the CsSnRK2 protein interacted with a sucrose non-fermenting protein (XP_004168589.1). We speculated that this protein belongs to a protein family other than SnRK2 and plays an important role. These results provide useful information for further analysis to verify the function of CsSnRK2 genes.