Characterization and expression analysis of the SPL gene family during floral development and abiotic stress in pecan (Carya illinoinensis)

Identification and characterization of SPL genes in pecan

The genome database, protein sequence, and pecan CDS sequences were downloaded from the Hardwood Genomics Project (https://www.hardwoodgenomics.org/) (Huang et al., 2019). The protein sequences of Arabidopsis, rice, and poplar SPL gene families were downloaded from the Plant Transcription Factor Database (PlantTFDB). The protein domain sequences of SPLs from the other three species were used as query sequences to blast the pecan proteome data with Local Blast P retrieval, and the E value was set to 1e−6. The blast results were sorted out to remove the repetitive sequence, and the NCBI-CDD (https://www.ncbi.nlm.nih.gov/cdd) and online software Pfam (http://pfam.xfam.org/) were used to further verify. The physical and chemical properties were analyzed using the online software Expasy (https://web.expasy.org/protparam/), and subcellular localization was predicted using Wolf psort (https://www.genscript.com/tools/wolf-psort).

Classification and phylogenetic tree construction of the SPL gene family

Bioedit software and the conserved domain sequences of SBP proteins in pecan were used to carry out multiple sequence alignments. For the evolutionary relationship, the candidate SPL protein domains of pecan, rice, Arabidopsis, and poplar were analyzed for multiple sequence alignment using ClustalW. Then, the comparison results were analyzed using MEGA 5.0 with a maximum likelihood (ML) (bootstrap=500) to construct the phylogenetic tree. Another phylogenetic tree was constructed using the SPL domains of all pecan protein sequences for further analysis via the same method.

Gene duplication and Ka/Ks value calculations

The collinearity (synteny) in the pecan, and between pecan and Arabidopsis, were analyzed using MCScanX software, and the E-value was set as 1  × 10⁻⁵, in accordance with previous studies (Wang et al., 2012). The ratio of nonsynonymous to synonymous nucleotide substitutions (Ka/Ks) was evaluated among segmentally duplicated gene pairs to detect the selection mode. First, the nucleotide coding sequences (CDS) of the CiSPL genes were aligned using the Muscle program with default options in MEGA 5.0 according to the protein sequence alignments. The Ka/Ks was then calculated based on multiple sequence alignments using MEGA 5.0 (Hu et al., 2019).

SPL gene structure and conserved motif analysis

GSDS software was employed to analyze gene structure. Simultaneously, we used the online software MEME (http://meme-suite.org/tools/meme) to predict and analyze the protein conservative motif of the SPL protein sequences. The MEME parameters were set as follows: maximum number of motifs, 20; motifs width: 6-60.

Category and number of cis-acting elements in CiSPL promoters

We extracted 1,500 bp promoter regions upstream of the translation start codon (ATG) of the CiSPL genes from the pecan genome website, which we labeled as promoter fragments. The cis-acting elements of each CiSPL were identified using the online program PlantCARE (Lescot et al., 2002).

Expression profiles of CiSPLs at various stages of flower and fruit development

In order to explore the role of CiSPLs in pecan growth and development, a comprehensive expression profile analysis was carried out. Normalized expression levels during female flower development were retrieved from our previous RNA-Seq data (Wang et al., 2019), and the CiSPL gene expression data in fruit development was retrieved from RNA transcriptome data (Bioproject ID PRJNA435846). The calculated fragments per kb of exon per million mapped reads (FPKM) were used to normalize the gene expression values (Huang et al., 2019). Next, the log2(FPKM) values of the CiSPL genes were used to generate heat maps with the software MEV4.9.

Plant materials and stress treatment

The gene expression patterns of 32 CiSPLs were detected in the leaves of pecan seedlings at the Institute of Botany, Jiangsu province and the Chinese Academy of Sciences. The pecan seedlings were cultivated from the seeds with the same size of ‘Pawnee’ and grown in fields under natural conditions. At the four-true-leaf stage, uniform and healthy pecan seedlings were transferred into a greenhouse. After being cultured in Hoagland’s solution for seven days, pecan seedlings were subjected to a variety of abiotic stresses. We added 10% polyethylene glycol (PEG) 6000 to the solution for drought treatment, and 3‰ (w/v) NaCl to the nutrient solution for salt treatment (Mo et al., 2020). Young leaves from the stress-treated plants were collected 0, 6, 12, and 24 h after treatment. All samples were immediately frozen in liquid nitrogen and stored at −80 °C until RNA extraction. The untreated pecan seedlings were used as the control groups.

RNA extraction and real time PCR analysis

Total RNA was isolated using the Universal Plant RNA Extraction Kit (BioTeke, Beijing, China), and reverse transcribed to cDNA with the PrimeScript™ RT reagent Kit (TaKaRa, Kyoto, Japan). The gene-specific primers of 32 CiSPLs were designed using CDS sequences and Primer 5.0 software, and all the primers are listed in Table S1. qRT-PCR was performed on ABI StepOnePlus (Thermo Fisher Scientific, Waltham, MA, USA) with SYBR Green Realtime PCR Master Mix (Thermo Fisher). Each reaction volume was 20 µL, and included 1 µL of diluted cDNA template, 0.4 µL of each gene forward and reverse primer, 10 µL of SYBR Mix, 8.2 µL of ddH2O, and 0.4 µL of ROX Reference Dye. The PCR reaction procedure was as follows: 95 °C for 30s, followed by 40 cycles at 95 °C for 5 s, 60 °C for 34 s, and the dissolution curve reaction procedure (95 °C for 15 s, 60 °C for 1 min, and 95 °C for 15 s). The pecan actin gene was used as the internal reference gene (Zhang et al., 2016). The relative expression level for all CiSPL genes was calculated using the 2^−ΔΔCt method (Livak & Schmittgen, 2001). ANOVA Duncan analysis was performed for significant difference analysis using SPSS 21.0 software.

Identification and annotation of SPLs in pecan

Referencing the SPL protein sequences of Arabidopsis, rice, and poplar, we searched all assumptive SPL genes of pecan using Local BLAST and BioEdit software. Through Pfam and NCBI domain analysis, we identified 32 candidate SPL genes (Table 1). Previously, Wu et al. (2019) found 30 CiSPL genes in the pecan genome. However, in this study, we found two CiSPL genes (CiSPL4c, CIL1118S0019; CiSPL9d, CIL1312S0023) in addition to the previously identified 30 CiSPL genes. The SBP box domains of the 32 CiSPL proteins were aligned and analyzed using BioEdit software (Table. S1). All CiSPLs had two zinc finger-like structures and NLS fragments except for CiSPL4c, CiSPL7b, and CiSPL9d, which indicated that the SBP domain was conserved in pecan. Further, the CiSPLs were named by orthology with Arabidopsis SPLs. The characterizations of CiSPL proteins were analyzed using the ProtParam tool. The lengths of the 32 CiSPL proteins ranged from 108 (CiSPL4c) to 1,071 (CiSPL16a), with molecular weights ranging from 12.184 kDa to 118.955 kDa, and theoretical isoelectric points from 5.25 to 9.36. This range of variations means that different CiSPL proteins may operate in different microenvironments. Most CiSPL genes were located in the nucleus, only two genes (CiSPL4c and CiSPL7b) were located in the chloroplast.

Cluster analysis and phylogenetic tree construction of SPL proteins in pecan

In order to explore the evolutionary relationships of SPL proteins in pecan and other species, we used the ML method on the SBP domains of 95 SPL proteins from pecan (32), Arabidopsis (17), rice (18), and poplar (28) to construct the phylogenetic tree (Fig. 1). In order to better understand the phylogenetic relationship of the SPL members, we also performed multiple alignment and phylogenetic tree analyses on the core SBP domains of all CiSPLs and AtSPLs. The 32 CiSPL genes were divided into eight subgroups (from I to VIII), described previously, based on phylogenetic analysis and SBP domain alignment (Fig. 2). The subgroups II, VI, and VII had the largest number of members (six). Each subgroup contained at least one SPL gene from the other three species (Arabidopsis, rice, and poplar). Subgroup I contained one SPL gene from Arabidopsis and rice, and two from pecan and poplar, suggesting the functional conservation of the SPL7 gene in plants. CiSPL genes showed high similarity with the poplar orthologs, and had the most similar number of genes in each subgroup. The SBP domain was relatively conservative in different species.

Motif and gene structure analysis of SPL genes in pecan

The CiSPL gene structure was analyzed using GSDS 2.0 (http://gsds.cbi.pku.edu.cn/) (Fig. 3). The number of introns in the 32 CiSPLs ranged from 0 to nine. Almost half of the 32 SPL genes (15 CiSPLs) had two introns. There were nine introns in seven CiSPLs (mainly members of group I and II), CiSPL2b and CiSPL8b had three introns each, CiSPL4c and CiSPL9d had no introns, and the remaining six CiSPL genes contained only one intron each. Most of the CiSPL genes within the same subgroups showed similar gene structures in terms of intron number and exon length. CiSPL3a and CiSPL3b were in group VI, and the length of their protein sequences was shorter than 200 amino acids, while in subgroup II, the SPL protein sequence was longer than 1,000 amino acids.

Other conserved domains besides the SBP domain also played roles in protein function. Based on the MEME software analysis of 32 CiSPL protein sequences, 20 motifs were identified and numbered from 1 to 20 (Fig. 4, Table S2). The number of motifs ranged from one to 17 in each SPL protein. All CiSPL proteins (except for CiSPL4c, CiSPL7b, and CiSPL9d) contained motif 1 (Zn-2), motif 2 (Zn-1), and motif 3 (NLS). CiSPL7b contained only motif 3, while CiSPL4c and CiSPL9d only had motif 2. Five CiSPL proteins (CiSPL3a, CiSPL3b, CiSPL4a, CiSPL4b, and CiSPL5) contained only motif 1, motif 2, and motif 3, while other CiSPL proteins had more motifs. CiSPL12a, CiSPL12b, CiSPL12c, and CiSPL12d had 17 motifs, and CiSPL16a and CiSPL16b had 11 motifs. Members of the same clade showed similar gene structures.

Gene duplication and collinearity analysis of CiSPL genes

Gene duplication is both the foundation of gene function diversification and the main driving force of species evolution. In order to explore the expansion and evolutionary mechanism of the CiSPL gene family, MCScanX software was used to analyze the SPL gene replication mode. Twelve genes were involved in nine segmental duplication events (Table 2). However, there were no tandem duplication events, which suggests that segmental duplications were the main mechanism of pecan SPL gene expansion. In order to further understand the evolutionary and functional relationship between the pecan and Arabidopsis SPL genes, we analyzed their collinearity relationships. The 16 collinear gene pairs, which included 13 CiSPLs and 11 AtSPLs, were as follows: CiSPL3a/AtSPL3, CiSPL4b/AtSPL5, CiSPL6a/AtSPL6, CiSPL6b/AtSPL6, CiSPL7a/AtSPL7, CiSPL7b/AtSPL7, CiSPL9a/AtSPL9, CiSPL9a/AtSPL15, CiSPL9b/AtSPL9, CiSPL9b/AtSPL15, CiSPL12b/AtSPL1, CiSPL13b/AtSPL13A, CiSPL3b/AtSPL3, CiSPL16a/AtSPL14, CiSPL16a/AtSPL16, and CiSPL16b/AtSPL16. This suggested that most SPLs had orthologs in Arabidopsis. For each duplicate SPL gene pair, the Ka/Ks ratio was calculated using MEGA 5.0 software to study divergence times. The Ka/Ks values in the pecan genome were distributed between 0.2–0.5, while the Ka/Ks ratios between the Arabidopsis and pecan genomes ranged from 0.1–0.4 (Fig. 5). This indicated that the SPL genes had a strong purifying selection between the pecan and Arabidopsis genomes.

Analysis of cis-acting elements in the promoter regions of CiSPL genes

We predicted and analyzed the cis-acting elements in the 1,500 bp promoter sequences of the CiSPL genes in order to understand the gene regulation mechanism. Except for the basic CAAT and TATA boxes, most of the other cis-acting elements mainly included stress response, hormone response, growth, and development, which were related to different life activities (Fig. 6). There were two cis-acting elements (GCN4_motif and AACA_motif) involved in endosperm expression and located in the promoter region of five and one CiSPL genes, respectively. The seed-specific regulation element (RY-element) was found only in the promoter region of CiSPL5. Nine meristem expression regulatory (NON-box and CAT-box) genes were identified in eight CiSPL promoters, and 11 zein metabolism regulation elements (O2 site) were also identified in eight CiSPL promoters. Additionally, circadian control (circadian) and flavonoid biosynthetic (MBSI) regulatory elements were located in the regions of three and one CiSPL promoters, respectively.

There were 46 MeJA-responsive elements (CGTCA motif and TGACG-motif), 19 auxin-responsive elements (TGA-element, TGA-box, and AuxRR-core), 10 gibberellin-responsive elements (GARE-motif, TATC-box, and P-box), 18 salicylic acid responsive elements (TCA-element), and 64 abscisic acid responsive elements (ABRE) located in the promoter regions of 15, 12, 9, 14, and 23 CiSPLs, respectively. There was the largest number of ABA-responsive elements, followed by MeJA-responsive elements, suggesting that CiSPL s may have a strong response to MeJA and ABA hormones, although this requires further experimental verification. The CiSPL transcription factor is an important adversity gene, and the number of hormone response elements in the promoters of CiSPLs varies widely. Among these, the CiSPL25 gene has the largest number of hormone response elements in its promoter region (up to 15). The promoter region of the CiSPL20 gene contains only one GA-related element, and no other types of hormone response elements. Different CiSPL s have different response patterns to different hormones.

There are five different types of stress-related elements in the promoter region of CiSPL genes. Among them, the most numerous are the regulatory elements related to anaerobic stress (ARE). Fourteen low temperature responsive elements (LTR) and 20 drought induction elements (MBS) were identified in 11 and 14 CiSPL genes, respectively. There are also stress-related response elements (TC_rich repeats).

Expression profiles of CiSPL genes during flower and fruit development

In order to understand the function of CiSPL genes, we analyzed the expression patterns of identified CiSPLs using the transcriptome data of female flower development (Table S3). Figure 7 shows that two-thirds of CiSPL genes were expressed during female flower development.CiSPL4c and CiSPL16b showed a high expression level during the whole flower development stage (FB1 to FL2). Four CiSPL genes (CiSPL8a, CiSPL8b, CiSPL9c, and CiSPL9d) were only expressed during the FL2 stage. The CiSPL4b gene was particularly expressed in the FL1 stage. CiSPLs play different roles in different stages of female flower development, and there is differentiation of gene function. We also analyzed the expression of CiSPL genes across three stages of pecan fruit development (Fig. 8, Table S4). Surprisingly, about half of the CiSPL genes were not expressed during fruit development. CiSPL12a, CiSPL12b, and CiSPL16a genes were highly expressed throughout fruit development. CiSPL3a and CiSPL12d were mainly induced and expressed during the early stages of fruit development. The specific functions of these CiSPL genes need to be further verified in future experiments.

Expression analysis of CiSPL genes under drought and NaCl stress treatment

In order to explore the response of CiSPLs to drought stress, the expression profiles of 32 CiSPLs were examined by qRT-PCR (Table S6). Figure 9 shows that nearly half of the CiSPLs were found to show a drought stress response. In general, shorter genes were expressed more quickly. Four of the shorter CiSPLs (CiSPL4a, CiSPL4b, CiSPL6d, CiSPL8a, and CiSPL13e) were induced rapidly and reached their peak at 6 h. CiSPL7a, -7b, -12a, -13f, and -16b expression reached their peak 12 h after drought stress. CiSPL6c was the most responsive gene, with an expression 100 times greater than that of control. Thirteen CiSPLs (CiSPL-3a, -4c, -5, -6b, -9a, -9b, 12c, -12d, -13a, -13b, -13c, -13d, and -16a) were not induced or down-regulated during the whole process. Generally, CiSPLs play an important role in drought stress response.

Under the NaCl treatment, nearly half of the 32 SPL genes were induced (Fig. 10, Table S7). However, six CiSPLs (CiSPL2b, -4b, -7a, 12a, -13f, and -16b) were down-regulated during the whole stress treatment process. CiSPL4c was significantly expressed at 6 h. The transcription level of six genes (CiSPL3a, -3b, -4a, -6a, -6c, and -9c) reached the maximum at 12 h. CiSPL8a and CiSPL13e were slightly induced at 12 h. CiSPL6d gene expression reached its peak at 24 h (nearly 20-fold over the control). CiSPL9d was induced at 12 h and 24 h, and the expression level was nearly the same.