Genome-wide identification and characterization of the MADS-box gene family in Salix suchowensis

Datasets and sequence retrieval

All the latest version files related to the Salix suchowensis genome sequence that were used for the identification of MADS-box genes were downloaded from the website of the Bioinformatics Laboratory of the Information College of Nanjing Forestry University (https://figshare.com/articles/Willow_gene_family/9878582/1). Arabidopsis genomic data and 89 MADS-box sequences were downloaded from The Arabidopsis Information Resource (TAIR, http://www.arabidopsis.org/index.jsp) with the accession numbers reported by Parenicová et al., and the MADS-box protein data for rice were obtained from the Rice Genome Annotation Project (RGAP, http://rice.plantbiology.msu.edu/index.shtml) (Kawahara et al., 2013; Parenicova et al., 2003).

Identification and distribution of MADS-box genes in willows

The method used to identify proteins corresponding to the willow MADS-box genes was similar to that used for other species (Duan et al., 2015; Tian et al., 2015; Wei et al., 2014). Fasta and Stockholm format files for the MADS-box domains were retrieved from the Pfam database (release 31.0, http://pfam.xfam.org/) with the accession number ‘PF00319’ (Finn et al., 2016). To obtain potential proteins, an alignment of MADS-box seed sequences in the Stockholm format was generated by a tool in the HMMER programs (hmmbuild) to build an HMM model, and then the model was used to search all willow proteins using another tool (hmmsearch) with the default parameters (Eddy, 1998). BLASTp ( E-value = 1⁻³) was used to align the Fasta profile downloaded from the PFAM website with all willow protein sequences (Willow.gene.pep) (Camacho et al., 2009). The potential willow MADS-box genes were obtained by taking the intersection of the above two results. To validate the confidence of these genes, we used the SMART programme (http://smart.embl-heidelberg.de/) to confirm whether a MADS-box domain was contained in each candidate MADS-box protein (Letunic, Doerks & Bork, 2015). Genes that did not contain an entire MADS domain were removed to identify eligible MADS-box gene family members. In addition, we used the ExPasy tool (https://www.expasy.org/tools/) to calculate the lengths, molecular weights, and isoelectric points of these putative MADS-box proteins. Finally, all identified MADS-box genes were mapped onto willow chromosomes with an in-house Perl script (http://bio.njfu.edu.cn/willow_chromosome/BuildGff3_Chr.pl). The distribution of each MADS-box gene on the willow chromosomes was plotted using the MapInspect software (https://github.com/quyanshu/Willow-gene-family/blob/master/BuildGff3_Chr.pl), and these genes were renamed based on their chromosomal distributions.

Multiple alignment and phylogenetic analysis of the willow MADS-box genes

The sequence logo of the identified willow MADS-box genes was generated using the web-based application WebLogo3 (http://weblogo.threeplusone.com) with the default parameters (Crooks et al., 2004). To obtain the conserved MADS-box domains of these willow MADS-box genes, we employed the online tool SMART and the PFAM database and used ClustalX (version 2.1) to perform multi-sequence alignment of the MADS-box domains obtained from SMART (Larkin et al., 2007). The online tool BoxShade (http://www.ch.embnet.org/software/BOX_form.html) was then used to colour the resulting alignment.

In general, all Willow MADS genes can be divided into two categories (M-type and MIKC-type) through the PlantTFDB website (http://planttfdb.cbi.pku.edu.cn/). However, to obtain a better subgroup classification of these genes, a multiple sequence alignment including willow (SsMADS) and Arabidopsis (AtMADS) MADS-box proteins was performed using Muscle, and a NJ tree was built with MEGA 7.0 based on this alignment (Edgar, 2004; Jin et al., 2014; Kumar, Stecher & Tamura, 2016). A NJ tree was then established for all Arabidopsis MADS-box proteins to check the reliability of this method (Duan et al., 2015). A phylogenetic tree was constructed using a similar method with the identified SsMADS domains and 66 rice MADS-box core domains (OsMADS). Additionally, a phylogenetic tree was built based on the identified SsMADS proteins.

Subsequently, to enable better comparison of MADS-box genes in Salicaceae, a phylogenetic tree was established for all SsMADS and Populus trichocarpa MADS-box genes. The method was consistent with that described above.

Finally, the orthologues of each SsMADS gene in A. thaliana, rice and Populus were determined based on the phylogenetic trees of the MADS-box domains or proteins and the BLASTP programme results (bi-direction, best hit, E-value = 1e⁻²⁰) (Chen et al., 2007).

Gene structure analysis of the willow MADS-box genes

The intron-exon structures of the willow MADS-box genes were contained in our own assembled protein annotation file. After annotation information for all SsMADS genes was extracted using a Perl language script, an intron-exon structure diagram was obtained from the online tool GSDS (Gene Structure Display server, http://gsds.cbi.pku.edu.cn/) (Hu et al., 2015).

Multi-sequence and BLASTp alignments (E-value = 1e−20) were performed to obtain the similarities between these SsMADS genes. To estimate gene duplication events in the SsMADS genes, the following metrics were set: (1) the proportion of regions used for alignment of the longer gene should exceed 65% and (2) the similarity of the aligned regions should exceed 65% (Bi et al., 2016).

To better reveal the structural features of the SsMADS proteins, the online tool Multiple Expectation Maximization for Motif Elicitation (MEME, http://meme-suite.org/) was used to predict conserved motifs in the encoded SsMADS proteins (Bailey et al., 2006). The parameters were set to a repeat motif site of any number, a maximum number of motifs of 15, and a width of each motif ranging from 6 to 60 residues. The web-based software 2ZIP (http://2zip.molgen.mpg.de/) was used to verify whether these SsMADS proteins contained the Leu zipper motif, and other important conserved motifs, including LXXLL and LXLXLX, were searched manually (Bornberg-Bauer, Rivals & Vingron, 1998).

Expression analysis of the willow MADS-box genes

To obtain more information regarding the roles of MADS-box genes in willows, RNA-Seq data from the sequenced genotype were used to quantify the expression levels of MADS-box genes in five tissues from S. suchowensis. The BWA programme was used to map back the S. suchowensis RNA-Seq reads from five tissues (roots, stems, leaves, buds and skins) onto the SsMADS gene sequences, and the number of mapped reads for each SsMADS gene in RPKM (reads per kilo base per million mapped reads) was calculated manually and standardized using Log₂ RPKM (Li & Durbin, 2009; Wagner, Kin & Lynch, 2012). A gene expression profile heat map was drawn with Bioconductor (pheatmap package) (Gentleman et al., 2004).

RNA isolation and Real-time quantitative RT-qPCR

Total RNA was isolated from five frozen willow tissues using an RNA kit (RNAprep Pure Plant Kit, Tiangen, Beijing, China), the specific procedures can be found in the manufacturer’s instructions. The quality and concentration of different RNA samples were determined by a NanoDrop 2000 c spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and 1.0 percent (w/v) agarose gel electrophoresis. cDNA was synthesized from 1,000 ng of total RNA in a 20 µL reaction volume using PrimeScript™ RT Master Mix (TaKaRa, Dalian, China) according to the manufacturer’s instructions. The resulting cDNA was then diluted three-fold and stored at −20 °C for the subsequent RT-qPCR assays.

For gene expression quantification, using the Oligo 7 algorithm (https://en.freedownloadmanager.org/Windows-PC/OLIGO.html) to design specific primers for each SsMADS gene, primer details are listed in the Table S1. The expression of 6 SsMADS genes was verified using RT-qPCR with a total volume of 10 µL per reaction (1 µL of cDNA, 1 µL of the forward and reverse primers, 5 µL of SYBR Green mix (TaKaRa, Dalian, China) and 3 µL ddH₂O) and performed on a StepOnePlus™ System (Applied Biosystems). The reactions were performed under the following conditions: 95 °C for 30 s and 50 cycles of 95 °C for 15 s and 60 °C for 1 min. The specificity of the amplicon for each primer pair was verified by melting curve analysis. All experiments were performed in three biological replicates; each replicate being measured in triplicate. Relative expression levels were calculated using the 2^−ΔΔCt method, with the OTU-like cysteine protease gene (OTU) of S. suchowensis as reference gene.

Identification and characterization of the MADS-box gene family in S. suchowensis

Sixty-four MADS-box genes were obtained using the HMMER toolkit to search the Hidden Markov Model of the MADS-box DNA-binding domain in the willow whole-genome protein sequence. The accuracy of the results was verified through BLASTP and HMMER mutual verification. Subsequently, the potential MADS-box genes were submitted to the SMART website for further verification. Four genes were removed due to lack of a MADS domain, and the remaining 60 probable MADS-box genes were selected as MADS-box superfamily members.

To better understand the MADS domain of S. suchowensis, a sequence logo and a multiple alignment with 60 SsMADS domains were generated. Amino acids 3, 23, 24, 27, 30, 31, and 34 were highly conserved, which confirmed conservation of the MADS domain (Fig. S1).

As shown in Fig. 1, the structures of the type I and type II SsMADS genes were quite different, and the type II SsMADS genes were more conserved than the type I genes. The MIKCc subgroup was the most conserved type, and several conserved motifs, including RQVT and RIEN, were concentrated at the N-terminus. The similarities between types I and II mainly occurred in the central region near the C-terminus. For example, differences in the N-terminal amino acids in Physcomitrella patens were reported to determine the differences between type I and type II MADS-box genes, whereas MIKCc and MIKC* are distinguished by the C-terminus (Henschel et al., 2002). In general, the type II MADS-box genes of S. suchowensis, particularly the MIKCc subgroup, were more conserved.

Detailed characteristics, including the classification, chromosomal distribution, homologous genes, and related physicochemical properties, of the SsMADS genes are listed in Table 1. As shown in Table 1, these protein sequences ranged from 80 amino acids (SsMADS34) to 894 amino acids (SsMADS40), with an average of 277 amino acids. Furthermore, the range of isoelectric points (PIs) also showed a large fluctuation, from 4.44 (SsMADS23) to 10.33 (SsMADS34), and the molecular weights (MWs) ranged from 9.20 kDa (SsMADS34) to 98.51 kDa (SsMADS40). These findings reflect the high complexity of willow MADS-box genes.

Chromosome distribution characteristics of the willow MADS-box genes

Fifty-six of the 60 SsMADS genes were distributed on 19 putative willow chromosomes, and these genes were renamed SsMADS1 to SsMADS56 based on their locations on the chromosomes. Only four SsMADS genes (willow_GLEAN_10001835, willow_GLEAN_10001302, willow_GLEAN_10001292, and willow_GLEAN_10000968) could not be mapped onto any chromosome, and these were renamed SsMADS57, SsMADS58, SsMADS59, and SsMADS60, respectively. As demonstrated in Fig. 2, chromosomes (Chr) 1 and 2 contained the largest number of SsMADS genes (six genes per chromosome), followed by Chr7, Chr8 and Chr9 (five genes per chromosome). Four SsMADS genes were found on Chr3 and Chr10, and three were found on Chr4, Chr6 and Chr16. Additionally, three chromosomes (Chr14, Chr15, and Chr17) contained two SsMADS genes, whereas only one SsMADS gene was found on Chr5, Chr11, Chr12, Chr13, Chr18 and Chr19. ChrN indicated that genes were not mapped on any chromosome.

The distribution of the MADS-box genes was not random; instead, an enrichment region showed a relatively high density on some chromosomes or chromosome fragments. Previous studies showed that a single chromosome region within 200 kb that contained two or more genes could be defined as a gene cluster (He et al., 2012; Holub, 2001). Genes that are used in large amounts are clustered in the genome to facilitate the rapid synthesis of large numbers of transcripts, which is important for predicting the potential function of co-expressed or clustered genes in angiosperms. According to the present study, a total of 21 SsMADS genes in willows were clustered into 11 clusters and distributed on nine chromosomes (Fig. 2). Two gene clusters were found on Chr1, including four SsMADS genes; one gene cluster each was distributed on Chr2, Chr3, Chr4, Chr7, Chr8, Chr9, Chr14 and Chr17. Three SsMADS genes were distributed in the gene cluster on Chr3, whereas no gene cluster was found on the other ten chromosomes.

Classification of MADS-box genes in willows

To better classify these SsMADS genes, a phylogenetic tree (NJ tree) was constructed using 88 AtMADS proteins from A. thaliana and the 60 SsMADS proteins identified in the present study. Based on the phylogenetic tree and structural features of the MADS-box proteins, all 60 SsMADS genes could be divided into two main groups (type I and type II) (Fig. 3). A total of 22 members were classified as type I (M-type), and these were further classified into Mα, Mβ and Mγ, with 11, seven and four members each, respectively. The remaining 38 members were categorized as type II (MIKC-type), which included 32 MIKCc-type and six MIKC*-type members. Furthermore, a similar classification was obtained with the NJ tree established for the 60 SsMADS domains and 66 rice MADS domains (Fig. S2). To better investigate the role of MADS-box genes in Salicaceae, we constructed a phylogenetic tree using 103 poplar and 60 willow MADS domains (Fig. S3). Based on the NJ tree described above, we found that most of the MADS-box genes from willows and poplars were clustered into sister pairs (40 SsMADS genes, accounting for 66.7% of all willow MADS-box genes, such as SsMADS32-PtMADS12 and SsMADS37-PtMADS89) because they originated from a common ancestor.

In addition, we compared the number of willow MADS-box genes with those of the ancient tree species Ginkgo biloba. The G. biloba MADS-box genes were predicted using the same method used to predict the willow MADS-box genes. The results revealed that G. biloba contained only 26 MADS-box genes, which was quite different from the number found in the willow genome. The number of MADS-box gene family members of gymnosperms, such as the Pinus taeda, an angiosperm variety, as well as monocotyledonous plants, such as Zea mays and Oryza sativa, and dicotyledons, such as Malus domestica and Glycine max, were also analyzed (Table 2). The gymnosperm genome was larger, but the number of this gene family was much smaller than that of the angiosperms.

Orthologues of SsMADS genes in Arabidopsis, rice and poplars

In this study, orthologous SsMADS genes in A. thaliana, rice and poplar were identified through a phylogenetic analysis combined with a BLAST-based method (bi-direction best hit). Finally, 35 pairs of orthologous genes from willow and A. thaliana, 35 pairs from willow and rice, and 57 pairs from willow and poplar were identified. The 22 type I SsMADS genes had 20 pairs of orthologous genes in poplar and five in A. thaliana, whereas rice contained no orthologues of the 22 type I SsMADS genes. The 38 type II SsMADS genes had 37, 30 and 35 pairs of orthologous genes in poplar, A. thaliana, and rice, respectively. In addition, 12 SsMADS genes were found to have identical domains in poplars (SsMADS9, SsMADS14, SsMADS17, SsMADS23, SsMADS24, SsMADS26, SsMADS43, SsMADS46, SsMADS50, SsMADS51, SsMADS53 and SsMADS58), and these accounted for 20% of the total number of genes. Among these 12 genes, 11 were MIKC-type, and only SsMADS23 was Mα; in addition, all 11 MIKC genes were found to have orthologous genes with high similarity in Arabidopsis and rice. For example, the similarity between SsMADS14 and OsMADS7/45 was 98.33%, the similarity between SsMADS14 and AGL2/AGL9 was 100%, the similarity between SsMADS43 and AGAMOUS was 98.31%, and the similarity between SsMADS50 and AGL2/AGL9 was 100%.

We also found that the vast majority of SsMADS genes that did not have orthologous genes in Arabidopsis also had no orthologous genes in rice.

Exon-intron structures of the SsMADS genes

To gain insights into the structural diversity of willow MADS-box genes, we analyzed the exon-intron organization of the coding sequences of each willow MADS-box gene. A striking bimodal distribution of introns was observed in the Arabidopsis, cucumber and apple MADS-box family genes; the MIKCc and MIKC*(Mδ) genes contained multiple introns, whereas the Mα, Mβ, and Mγ genes usually had either no or a single intron (Hu & Liu, 2012; Parenicova et al., 2003; Tian et al., 2015). We found a similar finding in willow. In Fig. 4, the SsMADS gene phylogenetic tree and the corresponding exon-intron structures are shown in the left and right panels, respectively. Among the 38 MIKC-type members, 34 (89%) members contained at least four introns, and the maximum of 13 introns was detected in SsMADS40. Correspondingly, among the 22 M-type genes, most of the members had no intron (77%) or a single intron, especially the Mγ-type SsMADS genes, and none of these four genes had any introns. Regardless, we found seven introns in SsMADS6 and eight introns in SsMADS8.

The following interesting phenomenon was also observed: the number of introns in the six MIKC*-type willow MADS-box genes was quite varied. Among these genes, SsMADS40 contained 13 introns, SsMAADS26 contained 10 introns, SsMADS31 contained nine introns, SsMADS28 contained four introns, and SsMADS34 and SsMADS56 contained only one intron each.

Gene duplication events and conserved motifs in willows

Two or more adjacent homologous genes located on a single chromosome are considered tandem duplication events (TDs), whereas homologous gene pairs between different chromosomes are defined as segmental duplication events (SDs) (Bi et al., 2016; Liu & Ekramoddoullah, 2009). In this study, we identified a total of 12 homologous gene pair (including 24 SsMADS genes) duplication events. Among them, 20 genes were MIKC-type genes (18 MIKCc and two MIKC*), and the remaining four genes were classified as Mα (Table S2). Besides, among the 12 homologous gene pairs, two appeared to have undergone TDs, and ten participated in SDs.

The conserved motifs of the 60 MADS-box proteins were predicted by the MEME programme to better analyze the sequence characteristics and structural differences among these genes. A total of 15 conservative motifs were predicted, and named from Motif 1 to Motif 15 (Fig. 5, Table S3).

Among these, Motif 1 and Motif 3 were widely present in all SsMADS genes. These two motifs were MADS domains, and Motif 1 was the most typical MADS domain. Motif 2 was a highly conserved K domain motif that is essential for protein interactions between MADS-box transcription factors and was present in all MIKC-type SsMADS genes except SsMADS44 and SsMADS56. Interestingly, the K-box domain was identified in SsMADS44 using the SMART programme but was not found using MEME because the two programmes used different algorithms. Further observation revealed that the K-box domain of SsMADS44 consisted of only 53 amino acids, whereas most K-box domains in willows were 92–93 amino acids in length; this shorter length might have been due to loss of a portion of the gene during evolution, which resulted in its distinctive features. Overall, SsMADS genes of the same subgroup had similar motifs, and we speculated that they might have similar functions. A total of six basic leucine zipper (bZIP) motifs were found in five SsMADS (SsMADS9, SsMADS16, SsMADS18, SsMADS19, and SsMADS46) using 2ZIP, and these motifs play important roles in the expression and regulation of higher plant genes. The activation domain LXXLL motif and the inhibitory domain LXLXLX motif were also found in willow MADS-box genes. In general, a large number of motifs with different structures and functions were found in the willow MADS-box gene family, indicating that the MADS-box genes play a variety of important roles in the gene regulatory network of willows.

Expression profiles of willow MADS-box genes in different tissues

The expression profile heat map of 60 SsMADS genes drawn using R is shown in Fig. 6. As illustrated in Fig. 6 and Table S4, most of the MADS-box genes were expressed at low levels or not expressed in these five tissues; this pattern was similar to the expression patterns of the MADS-box gene family in Medicago truncatula, in which seven of the genes, including SsMADS3, SsMADS12, and SsMADS18, were not expressed in the five tissues (Zhang et al., 2014). In contrast, 26 SsMADS genes were expressed in all tissues, and eight genes, including SsMADS9, SsMADS16, and SsMADS23, were highly expressed. SsMADS9 exhibited the highest expression level in four tissues (root, stem, leaf and bud) and showed high expression in bark. The gene belonging to the highly conserved MIKCc type, which can be considered the housekeeping gene of S. suchowensis, participates in various growth and development processes. SsMADS37 exhibited the highest expression in bark but quite low expression in the other four tissues. Additionally, seven of the eight genes with higher expression were of the MIKC-type; six of these were of the highly conserved MIKCc type, and the remaining gene was of the MIKC* type. Overall, the total RPKM value of the SsMADS genes was 287 in root and higher than 400 in the remaining four tissues. Therefore, the expression of the SsMADS genes in root was significantly lower than that in the stem, leaves, buds and bark. Thus, the MADS-box gene family plays a major role in willow morphogenesis.

RT-qPCR of six MIKC-type SsMADS genes were performed to validate their expression in RNA-seq. The results suggested that the expression levels of these genes were basically consistent with that of transcriptome sequencing data (Fig. 7). Furthermore, we found an interesting gene, SsMADS44, which was highly expressed in the stem but expressed at extremely low levels or not expressed (root) in the other four tissues by transcriptome data. However, it was also found to be highly expressed in Bud in RT-qPCR. The reason for these divergences might be the difference in sample growth (Meng et al., 2019).