Identification and expression analysis of SBP-Box-like (SPL) gene family disclose their contribution to abiotic stress and flower budding in pigeon pea (Cajanus cajan)

Plant material, treatments and growth conditions

For 2 weeks, Cajanus cajan (L.) Millsp. seedlings were grown in plastic pots (3–4 seedlings/pot) filled with peat moss, in a growth chamber at a humidity level of 60% and a temperature of 22°C at night and 25°C in the day. The plants were exposed to light for 16 h, followed by 8 h of darkness. Hoagland solution (10 mL/day) was applied to provide the necessary nutrients for plant growth. The seedlings were given 100 mM NaCl for the treatment of salt stress after 2 weeks of germination. The Hoagland solution was used as a preliminary treatment for the plants. The drought-stressed plants were deprived of water supply for 7 days whereas control plants were given Hoagland solution. For RNA extraction, the plants were divided into three replicates and each replicate contained three plants. The samples were frozen immediately using liquid nitrogen and then kept at −80°C until further experimentation.

Genome-wide identification and phylogenetic analysis of SPL family in C. cajan

A total of 11 members of the SPL gene family in A. thaliana (AtSPLs) were identified from the TAIR database (https://www.arabidopsis.org/). The AtSPL protein sequences were used as a query to identify the homologues in C. cajan (CcSPLs) using BLASTp (Protein BLAST: search protein databases using a protein query, nih.gov) (Swarbreck et al. 2007). The retrieved sequences were checked for redundancy using the SIAS tool by the Immunomedicine group (Immunomedicine Group: Tools >> SIAS, ucm.es) and annotated for further analysis using atom (https://atom.io/). The results were manually verified using NCBI and Uniprot (https://www.uniprot.org/) (Sherry et al. 2001; The UniProt Consortium 2015). After investigating the protein sequences retrieved from different species, phylogenetic analysis was performed (Pasquier et al. 2016). The maximum likelihood tree was constructed using IQ Tree (http://iqtree.cibiv.univie.ac.at).

Characteristics of SPL protein

The physiochemical characteristics of SPL gene family proteins were predicted using an online tool developed by the Swiss Institute of Bioinformatics called ProtParam (https://web.expasy.org/protparam/). Motifs of the SPL gene family were identified using MEME (Multiple Em for Motif Elicitation) (https://meme-suite.org/meme/tools/meme) where site distribution was zero or one sequence per occurrence and number of motifs was selected to 10. Motif logos were also designed by MEME (https://meme-suite.org/meme/tools/meme) (Bailey et al. 2015). The subcellular localisation of proteins was projected using an online-based soft berry tool called ProtComp 9.0 (ProtComp, softberry.com). SPL protein domains were predicated using InterPro (https://www.ebi.ac.uk/interpro) and were aligned using DNAMAN software (https://www.lynnon.com/; Wang 2015).

Gene sequence analysis of CcSPL

Exon and intron structures of CcSPL and AtSPL genes were predicated using GSDS Gene Structure Display Server Tool (ver. 2.0, http://gsds.gao-lab.org/) (Hu et al. 2015).

Analysis of putative promoter region of CcSPL genes

The 2000 bp promoter and UTR (untranslated region) were identified using the NCBI genomics tool (https://www.ncbi.nlm.nih.gov/nuccore) and the cis-elements in CcSPL gene promoter regions were identified using an online server-based website called PLANT CARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/; Lescot et al. 2002).

Expression profile from RNA-seq data

To study the protein expression, the transcriptomic data on abiotic stress (SRP104063) and flower budding (SRP262441) were retrieved from the NCBI SRA (Sequence Read Archive) database (https://www.ncbi.nlm.nih.gov/sra). All the available runs were mapped to C. cajan using Galaxy Bowtie2 (map reads against reference genome) (https://usegalaxy.org/). A tool of Galaxy called Cufflinks was used to calculate gene expression levels based on FPKM values (fragments per kilobase of transcript per million mapped reads, https://usegalaxy.org/). The expression profile of both abiotic stress and flower budding experiment was incorporated using the software TBTOOLS (https://bio.tools/tbtools; Chen et al. 2020).

RNA extraction, cDNA synthesis and qPCR

Thermo Nanodrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) was used to measure the amount of total RNA, extracted from C. cajan leaf samples using the Trizol technique. The Maxima H-minus First-Strand Synthesis kit was used to synthesise complementary DNA (cDNA) from 1 μg of RNA. To prepare it for future analysis, cDNA was kept at −20°C. A qPCR analysis was carried out using an iTaq Universal SYBR Green Super-Mix and a qPCR detection system (CFX96 TouchTM RT PCR Detection System, Bio_Rad Labs, Hercules, CA, USA). Gene-specific primers were designed using the internet tool ‘Oligo Calculator’ (http://mcb.berkeley.edu/labs/krantz/tools/oligocalc.html/, accessed on 13 July 2022), which were then confirmed by the NCBI-primer BLAST algorithm (https://www.ncbi.nlm.nih.gov/tools/primer-blast/, accessed on 13 July 2022). Each gene’s expression was analysed in triplicate, and the CcGAPDH gene served as the housekeeping gene (Ali et al. 2018).

CcSPL protein interaction analysis

To study the protein–protein interaction a web-based tool called string (https://string-db.org/) was used. All the sequences of CcSPL protein were pasted into the search box one by one and the results compiled manually using Microsoft PowerPoint into a single file (Bultan et al. 2017).

3D protein structure prediction and docking analysis

Protein docking analysis of the CcSPL proteins was undertaken using multiple tools. The first was to predict the three dimensional (3D) structures of the proteins. For most of the proteins of the family, 3D structures were predicted using an online server-based tool I-TASSER (Iterative Threading ASSEmbly Refinement) (https://zhanggroup.org/I-TASSER/). CcSPL3 was predicted using QUARK, because it contains less than 200 amino acids (https://zhanggroup.org/QUARK/). QUARK proposes developing an algorithm for ab initio protein structure prediction and protein–peptide folding, to construct the right protein 3D model using only the amino acid sequence (Waqas et al. 2021). To study the protein–ligand interaction, the three ligands, i.e. ABA, GA3 and IAA, were used and were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/). To perform molecular docking on the predicted SPL proteins, PyRx was employed. The best docking results (with the lowest values) were selected based on the binding affinity of the molecules, and the results were visualised using software called Chimera (Kim et al. 2019; Halim et al. 2021).

Ethical approval

This article does not contain any studies with humans or animals performed by any of the authors. We sought appropriate permission from the farm or field owner during specimens’ collection and experimentation. We confirm that during the collection and execution of the experiment, the authors have complied with the International Union for Conservation of Nature (IUCN) Statement on Research Involving Species at Risk of Extinction and the Convention on the Trade in Endangered Species of Wild Fauna and Flora. All methods were performed in accordance with the relevant guidelines/regulations/legislation.

Identification and phylogenetic analysis of CcSPL genes

A total of 11 SPL genes were identified in C. cajan (Table 1). The identified genes were manually curated for accession numbers, gene IDs, physiochemical properties of protein sequences, domains and subcellular localisation. The amino acid sequence length of CcSPLs ranged from 125 (CcSPL4) to 1000 (CcSPL15). All the proteins were located in the nucleus and the theoretical isoelectric points (pI) ranged from 6.04 (CcSPL1) to 9.45 (CcSPL15). The Grand Average of Hydropathicity (GRAVY) for each SPL protein was less than 0 (Liu et al. 2019). A maximum likelihood tree was constructed using IQ Tree (http://iqtree.cibiv.univie.ac.at) based on 79 SPL gene family sequences from C. cajan, A. thaliana, Cicer arietinum, Glycine max, Phaseolus vulgaris and Vigna unguiculata. All members of the tree were placed into nine groups (G1–G9: Fig. 1). Each group has at least one of the 11 members of the CcSPL protein (Nguyen et al. 2015). Two paralogous gene pairs were identified in G. max, including CcSPL1/GmSPL1 (G1) and CcSPL8/GmSPL8 (G5). For P. vulgaris, one paralogous pair CcSPL3/PvSPL3.1 was identified in G4. Similarly, a paralogous pair CcSPL13A/VuSPL13A was identified in G9 and CcSPL6/AtSPL6 in G8 for V. unguiculata and A. thaliana, respectively. In G9, CcSPL2.1/CcSPL2.2 was also identified as a paralogous pair from C. cajan. Gene pair CcSPL15/AhSPL15 in G6 was an out-paralog in the context of VuSPL15.1/15.2, GmSPL15, PvSPL15, and in paralogs for CaSPL15.

Fig. 1.

Unrooted phylogenetic tree of SPL protein from Cajanus cajan (Cc), Arabidopsis thaliana (At), Cicer arietinum (Ca), Glycine max (Gm), Phaseolus vulgaris (Pv), Vigna unguiculata (Vu) and Arachis hypogaea (Ah). The tree was generated using IQ Tree with the maximum likelihood method. Numbers above the branches show bootstrap value. Different colour circles represent different plants.

In G2, CcSPL14 was identified as orthologous to members including GmSPL14.1/14.2, PvSPL14, VuSPL14.1/14.2; and CcSPL7 in G3 was orthologous to GmSPL7/PvSPL7/VuSPL7. In G4, CcSPL4 was identified as orthologous to VuSPL3 and PvSPL4 (Yang et al. 2008).

A sequence alignment of these motif regions was performed for SPL proteins and represented as a sequence logo (Fig. 2). The Zn-1 (Cys-Cys-Cys-His) was conserved in all CcSPLs except CcSPL7, where the ‘His’ was replaced with another ‘Cys’. The second Zn-binding site (Zn-2) was conserved (Cys-Cys-His-Cys) among all the CcSPL proteins. Moreover, the Nuclear Localisation Signal (NLS) was highly conserved among all the c protein sequences (Ioannidi et al. 2016).

Fig. 2.

Sequence alignment to identify conserved motifs. Multiple sequence alignment of CcSPL proteins was done using t-coffee and the logo was designed using Weblogo (Crooks et al. 2004).

Gene structure analysis for exon/intron arrangement

Gene structures for exons and introns are shown in Fig. 3. The introns and exons of CcSPL are grouped according to phylogenetic analysis. The number of introns were not the same in all genes of the family as they ranged from 1 to 10. For instance, CcSPL1 contains 10, CcSPL7 contains 9, CcSPL14 contains 8, CcSPL2.2 contains 5, CcSPL2.1 and CcSPL13A contain 4, CcSPL6 and CcSPL8.1 contain 3, CcSPL15 contains 2, and CcSPL4 and CcSPL5 contain only one intron, respectively. Furthermore, 17 members of AtSPL and 11 members of CcSPL were analysed for the presence of conserved motifs using MEME software (ver. 5.4.1). Ten conserved motifs were detected and shown in Fig. 3. Motif-1 and Motif-4 were highly conserved in all of the proteins of the family. Motif-2 was also conserved in almost all of the genes except CcSL4. Motif-9 was also highly conserved and present in AtSPL1, AtSPL12, CcSPL1, CcSPL14, AtSPL14, AtSPL7, CcSPL7, CcSPL2.2, CcSPL2.1, AtSPL2, AtSPL10 and AtSPL11. Motif-3, 5, 6, 7 and 10 were conserved in AtSPL1, AtSPL12, CcSPL1, CcSPL14, AtSPL14 and AtSPL7. Motif-8 was conserved in CcSPL2.2, CcSPL2.1, AtSPL2, AtSPL10, AtSPL11, CcSPL15, CcSPL8, AtSPL8 and AtSPL15 (Fig. 3).

Fig. 3.

Exon–intron structures were characterised using GSDS 2.0, as shown in Fig. 3. Exons or CDS (coding sequence) are represented by yellow rectangles with rounded corners, introns by black lines, and untranslated portions by blue rectangles (UTRs). Intron stages are represented by 0, 1 and 2. Motifs are also displayed in this figure.

Protein–protein interactions analysis of CcSPL

Protein–protein interactions analysis was studied using a web-based tool String (Fig. 4a). The associations of SPL proteins with other proteins predicted through this tool are meant to be specifically significant (Takahashi et al. 2000). No interaction was observed for CcSPL2.1 and CcSPL2.2. CcSPL1 interacts with PEX4 (involved in peroxisome biogenesis) and AT2G34990 (involved in zinc ion binding). CcSPL14 interacts with PEX4 (involved in peroxisome biogenesis) and AGL20 (involved in transcription factor activity), which interacts with CcSPL3, CcSPL4, CcSPL8, CcSPL13A and CcSPL15.In contrast, AGL8 and CcSPL7 interact with COPT1, COPT2 and COPT6 (proteins involved in copper ion transmembrane transporter activity), FHY3 and HY5. Interestingly, CcSPL3, CcSPL4, CcSPL6, CcSPL8, CcSPL13A and CcSPL15 are working in a closer relationship whereas, all were interacting with the protein AGL8 (a probable transcription factor that helps in early floral meristem identity) in a complex network, which can be seen in the interaction diagram (Fig. 4).

Fig. 4.

(a) Using the online server-based technology called STRING, the protein–protein interactions of the SPL proteins with other members of the same family and other proteins. (b) Regulatory cis-elements were predicted for eight members of CcSPL family for 2000 bp promoter and UTR regions. The distribution of cis-elements related to various stress and hormone responses is depicted. Distinct shapes and colours are used to symbolise different cis-elements.

Cis-acting regulatory elements of CcSPL genes

Cis-acting regulatory elements are the main molecular switches for transcriptional regulatory activities of genes that affect a wide range of biological processes, i.e. hormone response, development and abiotic stress response (Yamaguchi-Shinozaki and Shinozaki 2005). Twenty-three types of hormones and stress-related cis-acting regulatory elements were predicted in promoter sequences analysis of CcSPL genes (Fig. 4b) (Supplementary file 1). Elements include ABRE (Abscisic Acid-Responsive Element), ARE (cis-acting regulatory element of anaerobic induction), BOX-4 (potentially involved in light responsiveness), CGTCA-MOTIF (which was involved in the hormonal MeJA-responsiveness), G-BOX (a trans-acting regulatory element with functions in light responsiveness), LTR (involved in low-temperature-responsive functions), 02-site (a trans-acting regulatory element involved in zein metabolism regulation), TGACG-motif (a trans-acting regulatory element involved in the MeJA-responsiveness), GCN4-motif (a trans-acting regulatory element related to endosperm expression), CAT-Box (a trans-acting regulatory element) and CAAT-box (cis-acting regulatory element related to meristem expression). All eleven members of the SPL gene family had the AT1 motif, GATA motif, TCT motif, AE box, GA motif, GT1 motif, TCCC motif, chs-CMA1a, LAMP element, box I and box II; each of these elements was a component of a light-responsive element (Mongkolsiriwatana et al. 2009).

Gene expression of CcSPLs in salt-tolerant and salt-susceptible pigeon pea genotypes under control and salt-treated roots

To gain insight into the functional properties of CcSPL family members, an RNA-seq analysis was performed using data from the NCBI SRA (Sequence Read Archive) database. A total of four runs were observed under bio project PRJNA382795. The current study involves the analysis of gene expression patterns of C. cajan root in two cultivars; ICP 7 (salt-tolerant pigeon pea genotype) and ICP 1071 (salt-susceptible pigeon pea genotype). A total of 39.25 gigabytes of experiment data was uploaded to SRA (SRP104063). The experimental data were used to evaluate the PKRM values using the Galaxy server-based tool and produced almost 61.04 gigabytes of data with more than 30 301 923 spots. The PKRM values were collected for all the CcSPL genes and the heatmap was generated using TB tools (Chen et al. 2020). In the heatmap, rows represent gene names and columns represent experiment values. The blue colour represents the downregulated genes, yellow represents no change of expression and red represents an upregulation (Fig. 5a). In the susceptible cultivar, CcSPL6 and CcSPL8 were downregulated whereas, CcSPL14 and CcSPL15 were upregulated under salt stress conditions. Similarly, in the salt-tolerant cultivar, all the CcSPLs were significantly upregulated in response to salt stress, except for CcSPL1, CcSPL2.2 and CcSPL4 (which were downregulated). It is important to notice that CcSPL14 and CcSPL15 were upregulated under salinity stress for both ICP 7 and ICP 1071 cultivars (Fig. 5a). In another SRA experiment, transcriptomes of two phases of unopened flowering buds of an isogenic C. cajan line ICPB 2043 were analysed. A total of eight runs were observed under bio project PRJNA634080. The expression of CcSPL genes can be divided into two distinct groups as highly expressed and weakly detected. The CcSPL1 and CcSPL8 were the most expressed transcripts at Bud Stage 1. Similarly, CcSPL8 and CcSPL15 were highly expressed at Bud Stage 2 as compared to other members of this family (Supplementary file 2) (Fig. 5b).

Fig. 5.

Gene expression analysis. (a) RNA-seq based expression profiles of CcSPL genes in response to salt stress in salt-susceptible and tolerant cultivars. (b) RNA-seq based expression profiles of CcSPL genes during two stages flower bud formation. The expression values [log2(FPKM)] are shown by the coloured bar. Red denotes high expression, whereas yellow denotes low expression; blue indicates no expression. Rectangles of various colours are used to denote the various CcSPL gene groupings. Table S1 lists the FPKM values of the CcSPL genes in various tissues. (c) SPL gene family relative qPCR test in response to abiotic (salt and drought) stresses. To produce an objective average value, the experiment was triplicated. In untreated plants, each gene had a default expression value of 1. Bars have been placed on each column to illustrate the standard error.

Gene expression analysis of CcSPL genes

To validate the expression of CcSPLs, real-time PCR-based expression quantification was performed for CcSPL genes in response to salt and drought stress (Hanafy et al. 2013). In response to drought stress, CcSPL1, 2.2 and 4 were significantly upregulated. On the other hand, CcSPL7, 14 and 15 were downregulated. The expression of other CcSPLs genes were not significantly modulated under drought stress conditions (Fig. 5c). In response to salt stress, CcSPL3, CcSPL4, CcSPL6, CcSPL13A and CcSPL15 were upregulated whereas, the expression of CcSPL7 was downregulated by salinity. It was observed that CcSPL4 is upregulated in response to both salt and drought stress conditions. When compared to control, CcSPL7 was shown to be downregulated under salt and drought stresses (Fig. 5c).

3D protein structure prediction and docking analysis

Protein 3D structures for all CcSPL proteins were predicted through I-TASER and the energy minimisation step was performed using Galaxy Mod-refiner. Since CcSPL3 and CcSPL14 protein expression were significantly modulated in response to environmental stimuli (Fig. 5), these were selected for structure prediction and docking analysis. Moreover, ABA and GA3 are important phytohormones for plant growth and stress response. Therefore, these two were selected as ligands in docking analysis. After the docking analysis, results were also collected for other CcSPL proteins and presented in the supplementary data (Supplementary files 3 and 4). The Fig. 6 displayed the 3D structure of CcSPL3 (Fig. 6a) and CcSPL14 (Fig. 6b) as well as the ligands attached to them. In Fig. 6a it can also be observed that ABA was attached to CcSPL3 at two Arg sites, two Glu sites, and at an Asn site with binding affinity of −4.4 kcal/mol, GA3 was attached at Arg, three Gln sites, Phe, Glu, Ile, Lys and Val sites with a binding affinity of −5.5 kcal/mol, and IAA was attached at Asp, Lys, Val, Glu and Arg with a binding affinity of −4.5 kcal/mol. Similarly, Fig. 6b showed ABA was attached to CcSPL14 at Ser, Ile, Arg, Leu, Trg, Ser and Gly sites with a binding affinity of −6.6 kcal/mol, GA3 was attached at two Arg sites, Ile, and multiple Ser sites with a binding affinity of −6.4 kcal/mol and IAA was attached at Phe, two Leu sites, Val, Tyr and Asp sites with a binding affinity of −5.5 kcal/mol.

Fig. 6.

(a) The 3D structure of CcSPL3 with distinct binding affinities of ABA, GA3 and IAA connected to it. (b) The ABA, GA3 and IAA bound to CcSPL14 in its 3D structure with various binding affinities.