2.1 Identification of SPL genes in four Ipomoea species
To obtain SPL genes in Ipomoea species, two methods, BLASTP and HMM, were used for identification, and two online databases, SMART and ScanProsite, were used for confirmation. As a result, a total of 29, 27, 26, 23 SPL genes were identified in I. batatas (Ib), I. trifida (Itf), I. triloba (Itb), and I. nil (In), respectively. For consistency, the Ipomoea SPL genes were renamed according to their chromosomal location order in each genome (Supplementary Table S1). The numbers of SPL genes and their percentage of the total genes in each species were displayed (Figure 1a). We found that I. nil genome has the least number of SPL genes compared to the other three species.
Subcellular localization analysis showed that most of the Ipomoea SPL proteins (94, 89.52%) were localized in the nucleus (Supplementary Table S1), suggesting their critical role in regulatory functions. Besides, the physical and chemical properties of SPL proteins were significantly different within species, but exhibited similar patterns among the four species (Figure 1b-d, Supplementary Table S1). To be more specific for all of the Ipomoea SPL proteins, the number of amino acids ranged from 103 (InSPL9) to 1141 (IbSPL11), the molecular weight (MWs) varied from 12.01 (InSPL9) to 126.51 (IbSPL11) kDa, the isoelectric points (pIs) ranged from 5.40 (ItfSPL14) to 10.55 (InSPL15).
2.2 Comparative phylogenetic analysis of Ipomoea SPL genes
To explore the evolutionary relationship among the Ipomoea SPL genes, a rooted neighbor-joining phylogenetic tree was constructed using 105 SPL proteins from four Ipomoea species and 51 SPL proteins from other five plant species (A. thaliana, S. lycopersicum, O. sativa and Chlamydomonas reinhardii) (Figure 2, Supplementary Table S2). Based on the classification of SPLs from A. thaliana, S. lycopersicum and O. sativa [7-9], Ipomoea SPL genes were classified into eight clades, which were numbered from I to VIII (Figures 1, 2). All clades had at least one SPL gene in each Ipomoea species, indicating the conservation of SPLs across Ipomoea genomes. However, the number of SPLs in certain clades were highly variable among Ipomoea species, suggesting the diversity of SPLs in genus Ipomoea. Clade I was the smallest subfamily, containing only one member for each Ipomoea species. Clade IV contained the highest number of SPLs (> 26%) in genus Ipomoea, and it could be further divided into three subclades (IV-a, IV-b, and IV-c). Members of IV-b and IV-c subclades only contain SPL genes from Ipomoea species, and no homologs of other species were found, which indicates that the Ipomoea SPL genes in these two subclades were generated in an ancient ancestor of the Ipomoea lineage. In addition, the phylogenetic analysis also indicated that most IbSPL genes were closer to ItfSPL genes than either ItbSPL or InSPL genes, supporting the fact that I. trifida is the most closely related diploid to hexaploid sweet potato [29]. When the number of SPL genes in Ipomoea species were compared with other species, we found that the number in Ipomoea species (the average number of SPL genes in 4 Ipomoea species was 26) was greatly increased by almost 1.6, 1.4 and 2 times than that in A. thaliana (16), O. sativa (19) and S. lycopersicum (13), respectively. The results indicated that Ipomoea SPL genes have undergone extensive expansion after the speciation of S. lycopersicum.
2.3 Gene and protein structure of the Ipomoea SPL family
To explore the structural diversity of the Ipomoea SPL genes, intron/exon structure analysis were carried out (Figure 3). Gene structure illustrations showed the high variation in the number of exons for Ipomoea SPL genes, ranging from 2 to 15 (Figure 3a, Supplementary Table S1). For instance, IbSPL20 possesses the largest number (15) of exons, whereas most of the Ipomoea SPL genes in Clade VI contained the least number (2) of exons. Moreover, most Ipomoea SPL genes in the same clade exhibited similar gene structures, even if they belong to different species. For example, SPL members in clade I and II contained the most exons, ranging from 10 to 15, while SPL members in the remaining clades possess 2 to 7 exons (except IbSPL19). Our results suggested that exon/intron gain or loss had occurred during the course of the Ipomoea SPL gene evolution, which may result in the functional divergence of SPL genes.
To investigate the features of Ipomoea SPL proteins, the conserved domains were analyzed by multiple sequence alignment. The results showed that all of the SPL members had the SBP domain, which was composed of two non-interleaved zinc finger-like structures (Zn-1/2) and one nuclear localization signal (NLS) motif (Figure 3b, Supplementary Figure S1). Based on the alignments, the Zn-2 motif is more conserved than the Zn-1 and NLS motif in all Ipomoea SPLs analyzed, which was also observed in Rosaceae [30] and Oryza species [31] (Supplementary Figure S1). For example, the Zn-2 motif in all Ipomoea SPLs was a Cys‐Cys‐His‐Cys (C2HC) type (except IbSPL23); and the Zn-1 motif showed different signatures, which was Cys‐Cys‐Cys‐Cys (C4) type in clade I and Cys‐Cys‐Cys‐His (C3H) type in the remaining clades. Moreover, other conserved domains were identified in specific clades. For instance, SPLs in clade I and II possess an DEXDc domain (Supplementary Figure S2), which were involved in ATP-dependent DNA unwinding [32]; SPLs in clade II contain Ankyrin repeats (Supplementary Figure S3), which were considered to be significant for mediating protein-protein interactions [33].
To gain a better understanding of Ipomoea SPL protein characteristics, the MEME software was used to explore the motif compositions (Figure 3d, Supplementary Figure S4). It was found that Ipomoea SPL proteins within the same clade showed similar motif compositions, while Ipomoea SPL proteins in different clades exhibited distinct differences in motif composition. In detail, two motifs (motif 2 and 3) existed in all Ipomoea SPL proteins, which were actually part of the SBP domain (Supplementary Figure S1); four motifs (motif 4, 5, 6 and 7) existed in clade I and II, and motif 4 was actually DEXDc domain (Supplementary Figure S2); motif 8 existed in clade II, which were make-up of Ankyrin repeats (Supplementary Figure S3); motif 10 presented in clade IV with unknown function. In summary, our results revealed that some motifs only presented in the specific clade, indicating that the SPL genes may be significantly different in function between different clades.
2.4 Gene duplication, orthology relationship and selective pressure of Ipomoea SPLs
To investigate gene duplication modes for Ipomoea SPL genes, MCScanX [34] was used to perform gene collinearity analysis within each Ipomoea species. All of the Ipomoea SPL genes were estimated to be duplicated (no singleton mode was found), and segmental (51, 48.57%) is the dominant mode compared to the other duplication modes: dispersed (39, 37.14%), tandem (13, 12.38%) and proximal (2, 1.90%) (Figure 4ab, Supplementary Table S1, Supplementary Figure S5). Our results indicated that segmental duplication has played a predominant role in the evolution and expansion of Ipomoea SPL genes. In addition, we found that tandem duplication was the most frequent event in IV-b and IV-c subclades, suggesting that SPL genes in these two clades were expanded by tandem duplication.
To assign the orthology relationships of the SPL genes, OrthoMCL [35] was used to detect orthologous groups across O. sativa, A. thaliana, S. lycopersicum and the four Ipomoea species. A total of 25 (1, 2, 2, 8, 3, 5, 2 and 2) orthologous groups in the eight clades (clade I to VIII) were identified, respectively (Supplementary Table S3). Among these groups, nine groups had genes start from O. sativa, suggesting that SPL genes in these groups may have occurred prior to the split of monocots and dicots; four groups had genes from A. thaliana but absent in O. sativa, implying that they originated after the divergence between monocots and dicots; ten groups had genes only existed in genus Ipomoea, indicating that they may originated in the common ancestor of Ipomoea lineage. Furthermore, the potential functions of some Ipomoea SPL genes could be inferred from their orthologs in O. sativa, A. thaliana and S. lycopersicum species.
To understand the divergence of Ipomoea SPL genes, the Ka, Ks and Ka/Ks ratio for all orthologous groups were calculated using PAML software [36] (Figure 4c, Supplementary Table S4). As a result, the mean Ka/Ks values for all clades were lower than 1.0, suggesting that Ipomoea SPL genes were evolved under the pressure of purifying selection. Genes in clade VIII showed the lowest mean Ka/Ks values (0.15) than those in other clades, indicating that they may evolve under the strong positive selection. In contrast, genes in subclade IV-b and IV-c exhibited the highest Ka/Ks values, implying that these two subclades have generally diverged much more rapidly than other clades.
2.5 miR156 target site of Ipomoea SPL genes
Using the publicly available miRNA transcriptomes (Supplementary Table S5), a total of 9 IbmiR156 members were identified in I. batatas (Figure 5a). To explore the roles of miR156-mediated post-transcriptional regulation of SPLs in the genus Ipomoea, the transcripts of all the 105 Ipomoea SPL genes were searched for the target site of miR156 using psRNATarget [37]. As a result, a total of 69 SPL genes was found to be the potential targets of miR156, including 10 InSPLs, 19 ItbSPLs, 18 ItfSPLs and 22 IbSPLs (Figure 5b, Supplementary Table S1). Among the miR156 target SPL genes, most of which (84 %) the sites recognized by miR156 were located downstream of SBP domain in the CDS region, then followed by the 3′-UTR (Figure 5b).
To further validate the predicted miR156-SPL interactions in I. batatas, the publicly available degradomes were used (Supplementary Table S5). The results showed that seven miR156-SPL interactions predicted by psRNATarget were confirmed by degradome sequencing (Figure 5c, Supplementary Table S1). To be more specific, the miR156 target sites of IbSPL9, IbSPL10, IbSPL12, and IbSPL15 were located in the CDS region, while the target sites of IbSPL16, IbSPL17, and IbSPL28 were located in the 3′-UTR region.
2.6 Cis-acting elements in promoters of Ipomoea SPL genes
To understand the regulatory mechanism and potential function of Ipomoea SPL genes, cis-acting elements were analyzed in the 2000 bp upstream sequence from the start codon for all SPL genes by using PlantCARE database [38]. A total of 4088 putative cis-acting elements were identified, and then were divided into four categories: light responsiveness, plant growth, phytohormone, and abiotic/biotic stress response (Supplementary Table S6, S7). As shown in Figure 6a, SPL genes in the same clades exhibited similar cis-acting element compositions in the promoter, indicating their conserved biological functions during the evolution in the Ipomoea. Among these categories, abiotic/biotic stress response category covered the largest portion (43.66%), followed by the light response (25.93%), phytohormone response (20.77%), and plant growth (9.64%) category (Figure 6b). In the abiotic/biotic stress response category, elements such as MYB/MYC (responds to abiotic stress signals), STRE (metal-responsive element), WUN-motif (wound-responsive element) and LTR (low-temperature-responsive element) were found. In the light responsiveness category, a series of elements were identified, such as Box 4, G-box, GT1-motif, TCT-motif, GATA-motif, and MRE, of which Box 4 motif was the most common (24%). As for the phytohormone response category, the ABRE element (responds to ABA), CGTCA-motif (responds to MeJA), ERE (responds to ethylene), TCA-element (responds to salicylic acid), and as-1 (responds to auxin) were the most common cis-acting elements, appeared in more than 60 Ipomeoa SPL genes. In the plant growth category, ARE element essential for the anaerobic induction, CAT-box related to meristem expression, and O2-site involved in zein metabolism regulation were the three major elements. In conclusion, the cis-acting elements analysis suggested that Ipomoea SPL genes participated in various biological processes.
2.7 Expression profiles of IbSPL genes in different tissues
Among the four Ipomoea species, I. batatas is the most important cultivated crop worldwide. To explore the putative roles of the IbSPL genes, the tissue-specific expression patterns of IbSPLs were analyzed in 8 tissues (four aboveground and four underground tissues) of two sweet potato cultivars (Xuzi3 and Yan252) using publicly available transcriptomic data (Supplementary Table S8) [39]. FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values were calculated and then used to evaluate gene expression levels (Supplementary Table S9). As shown in Figure 7a, the expression patterns of IbSPLs could be classified into three groups. The first group included 7 IbSPL genes (IbSPL6/IbSPL15/IbSPL22/IbSPL23/IbSPL24/IbSPL25/IbSPL26), with lowest expression levels (log2(FPKM) < 2) in all tissues. The second group was expressed at relatively high levels in all tissues, which was composed of 5 IbSPL genes (IbSPL4/IbSPL8/IbSPL11/IbSPL17/IbSPL20). The remaining 17 IbSPL genes belonging to the third group, which were highly expressed in some aboveground tissues, especially in shoot or young leaf. In addition, the gene expression profiles in underground tissues (fibrous and tuberous root) were investigated, and a total of 7 IbSPL genes were found to be expressed at a high level in underground tissues (mean FPKM > 10), such as IbSPL1, IbSPL4, IbSPL11, IbSPL17, IbSPL20, IbSPL21 and IbSPL28, implying their potential functions in root development. Furthermore, our results showed that IbSPL genes within the same clades exhibit distinctive expression patterns, such as IbSPL14, IbSPL17, IbSPL28, and IbSPL29 in clade VI.
To further investigate the expression profiles of IbSPL genes in underground tissues, publicly available transcriptomes of eight different stages during root development from cultivar “Beauregard” were used (Figure 7b, Supplementary Table S8, S9). The results showed that overall expression patterns of IbSPL genes in roots of cultivar “Beauregard” were similar to that in roots of cultivar “Xuzi3” and “Yan252”. To be specific, the expression levels of IbSPL17, IbSPL28, and IbSPL29 were the highest in storage root compared with undifferentiated and fibrous root; the expression level of IbSPL1 was the highest in undifferentiated root; the expression level of IbSPL10 was the highest in fibrous root; the expression levels of IbSPL4, IbSPL8, IbSPL11, IbSPL20, and IbSPL21 had no distinct variation in all tissues. Our results implied that these genes may play important roles in root development.
To confirm the expression patterns of IbSPLs derived from transcriptomic data, qRT-PCR analysis was performed for 11 selected IbSPL genes (IbSPL1, IbSPL4, IbSPL5, IbSPL12, IbSPL16, IbSPL17, IbSPL20, IbSPL21, IbSPL27, IbSPL28, and IbSPL29) in 11 tissues of cultivar “Xuyu34” (Figure 8). The results showed that the expression patterns of IbSPL genes were consistent between transcriptomic data and qRT-PCR results. We found the expression levels for IbSPLs were quite different in various tissues. For example, IbSPL27 and IbSPL29 were highest expressed in young leaf; IbSPL12, IbSPL16, and IbSPL21 were highest expressed in flower. Moreover, the expression levels of IbSPL17 and IbSPL28 were gradually increased during root development, indicating their differential roles in storage root development in sweet potato.
2.8 IbSPL genes in response to exogenous phytohormones
Phytohormones are known as regulators in plant growth and developmental processes, and the promoter analysis revealed that IbSPL genes might be regulated by various phytohormones. To reveal the potential roles of the IbSPLs in hormonal signaling pathways, qRT-PCR was used to analyze the expression of the 6 IbSPL genes (IbSPL1, IbSPL4, IbSPL16, IbSPL17, IbSPL21 and IbSPL28) under exogenous phytohormone treatments, including indole-3-acetic acid (IAA), methyl-Jasmonate (MeJA), zeatin (ZT), and abscisic acid (ABA). In general, expression analysis indicated that the 6 IbSPL genes exhibited highly divergent response patterns under the phytohormone treatment in root (Figure 9). Under IAA treatment, IbSPL16 and IbSPL21 were rapidly upregulated after 6 h of IAA treatment, while the other four genes were upregulated after 12 or 24 h of IAA treatment. Under MeJA treatment, all IbSPL genes were significantly upregulated after 48 h of MeJA treatment, of which IbSPL4 and IbSPL21 were upregulated 45.3 and 178.1 folds. Under ZT treatment, all IbSPL genes were highly upregulated after 6 h of ZT treatment, then gradually declined, and finally were expressed at least 10 folds than CK. In particular, IbSPL21 showed extremely positive response to ZT treatment. Under ABA treatment, all examined IbSPL genes, particularly IbSPL1, IbSPL4, IbSPL16, IbSPL17, and IbSPL21, were rapidly upregulated after 6 h of ABA treatment, whereas IbSPL28 was significantly upregulated after 24 h of ABA treatment.
2.9 Regulatory sub-networks involving IbSPLs and other I. batatas genes in the storage root of sweet potato
To identify the regulatory sub-networks involving IbSPLs in the storage root of sweet potato, weighted gene co-expression network analysis (WGCNA) was performed based on transcriptomic data of mature storage root from 88 sweet potato accessions (Supplementary Figure S6, Supplementary Table S11). A total of 19 modules were obtained, among which three modules (turquoise, blue and yellow) were found contain 8 IbSPL genes in total (Figure 10a). In the yellow module, IbSPL20 had only one co-expressed gene. In the blue module, IbSPL9 and IbSPL10 were co-expressed with 17 and 107 genes, respectively. In the turquoise module, IbSPL1, IbSPL16, IbSPL17, IbSPL21, and IbSPL28 had 12, 147, 370, 732, and 198 co-expressed genes, respectively (Supplementary Table S12). Interestingly, IbSPL16, IbSPL17, IbSPL21, and IbSPL28 in the turquoise module share 101 co-expressed genes (Figure 10b), indicating they may belong to the same biological process.
To further explore putative functions of the SPLs in the storage root, GO enrichment analysis was performed for the co-expressed genes. For the IbSPL1, IbSPL20 and IbSPL9, we can’t obtain GO enrichment results due to the small number of co-expressed genes. For the IbSPL10 in the blue module, co-expressed genes related to “response to chitin” and “response to organonitrogen compound” were enriched (Supplementary Table S13). For the IbSPL16, IbSPL17, IbSPL21, and IbSPL28 in the turquoise module, the similar GO terms were enriched, such as “regulation of root morphogenesis”, “cell division”, “cytoskeleton organization”, “plant-type cell wall organization or biogenesis”, and “cellulose biosynthetic process” etc. (Figure 10c, Supplementary Table S13). It is known that the cytoskeleton, cell division, and cell wall organization/biogenesis are important biological processes in the storage root formation and development [40]. Our results indicated that IbSPL16/IbSPL17/IbSPL21/IbSPL28 may play a key role in storage root development in sweet potato.