Molecular phylogeny of mulberries reconstructed from ITS and two cpDNA sequences

Plant materials

We selected 33 mulberry accessions (seven species: M. alba, M. mongolica, M. cathayana, M. wittiorum, M. nigra, M. yunnanensis, and M. notabilis) for this study (see Table 1). Morus yunnanensis was obtained from the Institute of Sericulture and Apiculture, Yunnan Academy of Agricultural Sciences, Mengzi, Yunnan province, China. The other mulberry accessions were obtained from the Mulberry Germplasm Nursery at Southwest University, China (He et al., 2013; Zeng et al., 2015), where they are preserved by propagation through grafting. The relatively closely related species Artocarpus heterophyllus was selected as the outgroup for the phylogenetic analyses based on ITS sequence (KT002551) and cpDNA sequences (MG434693).

DNA preparation, sequence selection, and amplification

Genomic DNA was extracted from all mulberry accessions using the CTAB method (Saghai-Maroof et al., 1984). Chloroplast DNA was extracted as described by Shi et al. (2012). Based on the alignment of six whole chloroplast genomes (from M. notabilis (KP939360.1), M. indica (DQ226511.1), M. cathayana (KU981118.1), M. mongolica (KM491711.2), M. alba var. multicaulis (KM491711.2), and M. alba var. atropurpurea (KU355276.1)) using VISTA viewer (http://genome.lbl.gov/vista/index.shtml) (Frazer et al., 2004), we selected the two variable cpDNA sequences, trnL-trnF and trnT-trnL, for analyses. The ITS and the trnL-trnF sequences were amplified using primers described elsewhere (Taberlet et al., 1991; Zeng et al., 2015). We designed primers to amplify trnT-trnL (F: 5′-TGCGATGCTCTAACCTCT-3′; R: 5′-TAGCGTCTACCAATTTCG-3′). These sequences were amplified using GoTaq Flexi DNA Polymerase (Promega Corporation, Madison, WI, USA) according to the manufacturer’s instructions. All the polymerase chain reaction (PCR) cycles consisted of initial denaturation of 95 °C for 5 min; followed by 32 cycles of 95 °C for 30 s, annealing for 30 s, 72 °C for 1 min; and then final extension for 7 min. The amplified fragments were isolated by electrophoresis on 1% (w/v) agarose gels. The purified sequences were cloned into the pMD19-T vector, and then 7–20 positive clones were sequenced for each accession.

Sequence analyses

The cpDNA and ITS sequences were assembled and corrected for PCR errors using Sequencher 4.2 (Gene Codes Corp., Ann Arbor, MI, USA). All the sequences were aligned using Clustal X 1.81 software (Thompson et al., 1997). The boundaries of the ITS1, 5.8S, and ITS2 regions were determined as described elsewhere (Zeng et al., 2015). The GC content and sequence length was calculated using BioEdit (Hall, 1999). Three angiosperm conserved motifs (motif 1: GAATTGCAGAATCC, motif 2: TTTGAACGCA, motif 3: CGATGAAGAACGTAGC) were detected by BioEdit (Kerbs et al., 2017; Yakimowski & Rieseberg, 2014). The minimum free energy of the secondary structure was predicted using Mfold (http://www.bioinfo.rpi.edu/applications/mfold) (Zuker, 2003).

Phylogenetic analyses

Phylogenetic analyses of ITS sequences and cpDNA sequences were conducted using maximum-likelihood (M-L) and Bayesian inference (BI) methods. The first phylogenetic analysis was conducted using the ITS sequences of the 33 mulberry accessions determined in this study and those from another nine species (M. mesozygia (HM747171), M. insignis (HM747169), M. serrata (HM747176), M. rubra (HQ144180), M. celtidifolia (HM747168), M. macroura (HM747170), M. mongolica (KF784879), M. wittiorum (AY345154), and M. australis (KT002555)) reported previously (Zeng et al., 2015). The best-fit model SYM+G for BI analyses of ITS sequences was selected by the lowest Akaike Information Criterion (AIC) scores in MrModelTest 2.3 (Nylander et al., 2004). MrBayes v3.2.6 software was used for BI analyses (Huelsenbeck & Ronquist, 2001). Four Markov chain Monte Carlo chains were run for 2,000,000 generations, with sampling every 100 generations. The first 5,000 trees were discarded as burn-ins and the 50% majority-rule consensus tree was determined to calculate the posterior probabilities for each node. The standard deviations of split frequencies were checked and the number of minimum generations required for analysis were those with a standard deviation value lower than 0.01. The M-L phylogenetic trees were constructed using the most suitable Kimura 2-parameter model with the lowest Bayesian Information Criterion scores in MEGA 7 (Kumar, Stecher & Tamura, 2016). A discrete gamma distribution was used to model evolutionary rate differences among sites (five categories (+G, parameter = 1.3541)) (Kumar, Stecher & Tamura, 2016). All characters were equally weighted.

Phylogenetic analyses of cpDNA sequences were performed using sequence matrix data for the trnL-trnF and trnT-trnL regions. The best-fit model (GTR) for BI and M-L analyses for cpDNA sequences was selected by the lowest AIC scores in MrmodelTest 2.3 (Nylander et al., 2004). The other parameters for BI analyses based on cpDNA sequences were the same as the parameters for analyses based on ITS sequences. The M-L phylogenetic tree was constructed in MEGA 7 with 1,000 bootstrap replicates (Kumar, Stecher & Tamura, 2016). All characters were equally weighted.

Variations among ITS sequences

Thirty-three mulberry accessions were sequenced and their ITS sequences were analyzed (Table 1). For each accession, 7–20 clones were sequenced (Table 1). Of the 33 mulberry accessions, 19 contained only one type of ITS sequence, with lengths ranging from 611 bp to 631 bp. The ITS sequences of M. notabilis, M. yunnanensis, and M. nigra were confirmed to have lengths of 631 bp, 631 bp, and 624 bp, respectively, as reported previously (Zeng et al., 2015). Of the 19 mulberry accessions with one ITS sequence, 16 had an ITS sequence of 611 bp, with differences at only two single nucleotide polymorphism (SNP) sites.

The other 14 mulberry accessions had polymorphic ITS sequences (Table 1). There were two types of ITS sequences: short ITS sequences (ITS-α) and long ITS sequences (ITS-β). The ITS sequences of the mulberry accessions Shimiansang, Hanguodabaizhenzhu, Shuisang, and Pisang2 contained two short-type sequences. Sequences with higher and lower copy numbers were designated as ITS-α1 and ITS-α2, respectively. The accession Yun6muben contained two long ITS sequences (625 bp and 624 bp, designated as ITS-β1 and ITS-β2, respectively). The ITS-α sequences in 13 mulberry accessions were identical with a length of 611 bp. The length of ITS-β sequences ranged from 624 bp to 626 bp. An alignment of the ITS sequences of the mulberry accessions is shown in Fig. 1. We detected the 13-bp InDel reported previously (Zeng et al., 2015). Further, we found three types of 13-bp sequences among the 33 mulberry accessions (CGTATACAATGCG, TGTGTGCAATGCG, and CGTACACAATGCG). Alignment analyses of the ITS sequences revealed six 1-bp and three 2-bp InDels. Other sequence variations were SNPs.

Identification of ITS pseudogenes

The GC content of the ITS region, the minimum free energy of the secondary structure, and conserved motifs in the 5.8S rDNA region were used to identify ITS pseudogenes. As shown in Fig. 2 and Table S1, the GC contents of Shuisang-α2, Jianpuzhai-β, and Yun6-β were similar to that of ITS-α1/α. The GC content of Yun6muben-β2 was lower than that of ITS-β1. The GC contents of ITS-α2/β2 sequences in the remaining mulberry accessions were lower than those of ITS-α1/β1 sequences in at least one region of ITS1, 5.8S rDNA, and ITS2. The minimum free energy of secondary structure showed a similar pattern to that of GC contents, except that it was lower in Shuisang-α2, Jianpuzhai-β, Yun6-β, and Yun6muben-β2 than in ITS-α1/α and Yun6muben-β1. The sequences of Hanguodabaizhenzhu-α2, Shuisang-α2, Yun7-β, Jianpuzhai-β, Yun6-β, and Yun6muben-β2 contained all three conserved motifs. Based on the sequence length information, the GC content of the ITS region, the minimum free energy of secondary structure, and conserved motifs in the 5.8S rDNA region, we identified 11 ITS sequences as pseudogenes (Yun6muben-β2, Pisang2-α2, Shimiansang-α2, Shanxitiansang-β, Xinyizhilai-β, Lunjiao109-β, Huosang-β, Agentingsang-β, Hanguodabaizhenzhu-α2, Baojing7-β, and Yun7-β).

ITS phylogenetic analyses

We included M. mesozygia (clade D1), M. insignis (clade D2), M. serrata (clade A4), M. rubra (clade A3), M. celtidifolia (clade A1), M. macroura (clade B), M. mongolica (clade B), M. wittiorum (clade B), and M. australis (clade B) in the phylogenetic analyses, as they represented six of the eight clades classified in a previous study (Zeng et al., 2015). We constructed M-L and BI trees. The main clades in the M-L trees (Fig. 3) were the same as those in the BI trees (Fig. S2), and both were consistent with the results of a previous study (Zeng et al., 2015). However, the pseudogene Yun6muben-β2 clustered together with the functional ITS sequences in the BI tree (Fig. S2), so we used the M-L tree for further analyses. Most of the mulberry accessions were clustered in clade B in these two trees. In clade B, all the functional ITS sequences were grouped into three clades (clade 1, clade 2, and clade 3) (Fig. 3): two short-type ITS sequence clades (clade 1 and clade 2) and a long-type ITS sequence clade (clade 3). The ITS pseudogenes were grouped into two pseudogene clades, designated as pseudogroup 1 and pseudogroup 2 (Fig. 3).

Characteristics and phylogenetic analyses of cpDNA sequences

Based on the synteny analyses of the whole chloroplast genomes of M. indica (Ravi et al., 2007), M. notabilis (Chen et al., 2016a), M. mongolica (Kong & Yang, 2016), M. cathayana (Kong & Yang, 2017), M. alba var. multicaulis (Kong & Yang, 2017), and M. alba var. atropurpurea (Hu et al., 2014) (Fig. S1) and sequence alignment analyses, the trnL-trnF and trnT-trnL regions were selected for the cpDNA phylogenetic analyses. In the 33 mulberry accessions, the length of trnL-trnF ranged from 913 bp to 931 bp, and that of trnT-trnL ranged from 1108 bp to 1150 bp. The combined data matrix of trnL-trnF and trnT-trnL consisted of 2109 aligned nucleotides (Fig. 4).

We constructed M-L and BI phylogenetic trees using the trnL-trnF and trnT-trnL data (Fig. 5; Fig. S3). The topologies of the BI trees were disordered (Fig. S3), and some identical sequences were clustered into different clades. Therefore, the M-L trees were used for further analyses. First, M. notabilis and M. yunnanensis were diverged first among the 33 mulberry accessions. The remaining 31 mulberry accessions were separated into two main clades. Clade I contained 11 accessions of species M. alba (Banqiao6, Basailuona, Huasang, Huai302, Huosang, Shanxitiansang, Shimiansang, Sililanka, Xinyizhilai, Yidachimu, and Wupisang), two accessions of species M.mongolica (Jimengsang and Mengsang) (Figs. 4 and 5; Table 1). Clade II contained 10 accessions of species M. alba (Agentingsang, Gailiang10, Hanguodabaizhenzhu, Jianpuzhai, Leshandahongpi, Lunjiao109, Shuisang, Taiwanchaochangguo, Xinjiaposijiguosang, and Zhenzhubai), three accessions of species M. cathayana (Baojing7, Gui23, and Pisang2), four accessions of species M. wittiorum (Yun6, Yun6muben, Yun7, and Yun7muben), and one accession of species M. nigra (Fig. 5; Table 1). There were more sequence variations in clade II than in clade I (Fig. 4). Species of M. alba (Hanguodabaizhenzhu, Lunjiao109, Shuisang, Taiwanchaochangguo, and Xinjiaposijiguosang) formed a subclade, and M. nigra was separated from the other mulberry accessions on a long branch (Fig. 5).

Characterization and mutation pattern of ITS pseudogenes

ITS pseudogenes have been detected in many plants (Fan et al., 2014; Won & Renner, 2005; Xiao, Moller & Zhu, 2010; Xu et al., 2017), and should not be excluded from phylogenetic analyses of ITS sequences (Bailey et al., 2003). However, ITS pseudogenes have evolved at a faster rate than functional ITS sequences, which can cause confusion in phylogenetic analyses (Bailey et al., 2003; Fan et al., 2014). Thus, comprehensive analyses of ITS sequences are required. The ITS pseudogenes can be easily identified based on their GC content, minimum free energy of secondary structure, presence of conserved motifs, substitution rates, phylogenetic positions, and copy numbers (Bailey et al., 2003; Queiroz, Batista & De Oliveira, 2011). In the present study, we found ITS pseudogenes in 11 mulberry accessions, implying that incomplete concerted evolution is occurring in mulberry. Compared with functional ITS sequences, the ITS pseudogenes of mulberry showed higher diversity in the whole ITS region (ITS1, 5.8S rDNA, and ITS2), as found in other studies (Fan et al., 2014; Xiao, Moller & Zhu, 2010). Both the pseudogenes and functional ITS sequences were a mixture of short- and long-type sequences. Recombination was detected by the RDP4 program (Martin et al., 2015). The results showed that the short-type pseudogene sequences Shimiansang-α2 and Hanguodabaizhenzhu-α2 have recombined with putative long-type ITS sequences, representing one pathway of concerted evolution of ITS sequences.

Phylogenetic analyses and evolution of ITS

We conducted phylogenetic analyses of ITS functional sequences and pseudogenes. The main topologies of the phylogenetic trees were the same as that reported elsewhere (Zeng et al., 2015). In the present study, ITS pseudogenes were identified and separated from the functional ITS clades (Fig. 3) (Bailey et al., 2003). Clade 3 is a newly identified clade with long-type ITS sequences compared with the previous study (Zeng et al., 2015). Based on the sequence similarities and sequence length of the ITS region, the long-type ITS sequences in clade 3 may have closer relationships with the sequences in clades A, C, and D (Fig. 1).

In our analyses, most of the sequences were clustered in clade 1 and clade 2, with only three sequences clustered into clade 3 (Fig. 3). We detected inconsistencies between traditional and molecular classifications in this study. Twenty-nine mulberry accessions with morphological polymorphisms clustered together with a low phylogenetic resolution and had similar or identical ITS sequences. For example, Yun6muben and Yun7muben showed similar fruit, leaf, bud, and internode morphologies, indicating that both should be classified as M. wittiorum. However, Yun6muben had long-type ITS sequences that clustered in clade 3, while Yun7muben had short-type ITS sequences that clustered in clade 2. As the morphological characteristics were very complex, no details are discussed here. These phenomena are indicative of concerted evolution or introgression of the ITS region, and of hybridization (Bailey et al., 2003; Xu et al., 2017). Thus, the classification of mulberry is very complex, and there is considerable confusion at present. Systematic classification of mulberry resources should be based on molecular markers, morphology, chromosome number, genome data, and other characteristics.

Intra-individual polymorphisms in the ITS region have been detected in many plants and pseudogene sequences have been included in phylogenetic analyses (Mayol & Rossello, 2001; Xiao, Moller & Zhu, 2010; Zheng et al., 2008). Bidirectional, unidirectional, or different rates of evolution have been detected from analyses of ITS sequences (Wendel, Schnabel & Seelanan, 1995; Xu et al., 2017; Zheng et al., 2008). In this study, we detected multiple functional ITS copies, putative ITS pseudogenes, or recombinants of ITS pseudogenes in the same individual. The successive evolution patterns of ITS sequences are indicative of concerted evolution (Figs. 1 and 3). There are several lines of evidence for the concerted evolution of ITS sequences. First, clades A, C, D, and 3 had long-type and ancient ITS sequences. Second, Yun6 was identified as the progeny of Yun6muben through natural pollination. Yun6-β showed one SNP variation compared with Yun6muben-β1, suggesting that concerted evolution started in the F1 generation. Third, Hanguodabaizhenzhu-α2 and Shimiansang-α2 were found to be recombined from long-type ITS sequences, and other ITS pseudogenes showed more variations. Fourth, Shuisang-α2 had only one SNP variation compared with ITS sequences in clade 1 and clade 2, implying that concerted evolution is continuing or has completed only recently. Finally, most functional ITS copies were short-type sequences. In conclusion, the evolutionary process of ITS, as indicated by our data, is that long-type ITS sequences are undergoing concerted evolution to form short-type ITS sequences.

Hybridization is an important evolutionary mechanism in plants, especially in flowering plants (Kerbs et al., 2017; Soltis & Soltis, 2009; Yakimowski & Rieseberg, 2014). Hybridization can also help to explain the concerted evolution of ITS sequences (Xu et al., 2017). We detected three hybridized mulberry accessions (Jianpuzhai, Shuisang, and Yun6), which contained two functional ITS sequences. Hybridization between M. rubra (A clade) and M. alba (B clade) has been detected in previous studies (Burgess & Husband, 2004; Burgess, Morgan & Husband, 2008). It has been suggested that hybridization between M. rubra and M. alba caused the local decline of M. rubra (Burgess et al., 2005). This may indicate that M. alba has higher ecological potential. The continuous back-cross hybridization of M. rubra with M. alba may have led to the disappearance of ITS sequences from M. rubra, or the concerted evolution to form short-type sequences. Like other plants, local Morus species face the risk of extinction (Burgess & Husband, 2004; Burgess, Morgan & Husband, 2008; Ellstrand & Schierenbeck, 2000; Wolf, Takebayashi & Rieseberg, 2001).

Polyploidy is another important evolutionary mechanism that is known to affect the concerted evolution of ITS sequences in plants (Ainouche et al., 2004). The ITS sequences can retain their subgenomic sequences or become homogenized through concerted evolution (Wendel, Schnabel & Seelanan, 1995). Existing mulberry species show various polyploidy levels, with 14, 28, 35, 42, 49, 56, 84, 112, 126, or 308 chromosomes (Xuan et al., 2017; Zeng et al., 2015). Multiple nucleolus organizing region (NOR) loci have also contributed to the concerted evolution of ITS sequences. Two pairs of NOR loci have been reported for M. notabilis (Xuan et al., 2017). The other mulberry accessions contained at least two pairs of NOR loci (data not shown). Mulberry is often cultivated by grafting or propagated from cuttings, resulting in a long generation time. This may be another factor affecting the concerted evolution of ITS sequences. In summary, hybridization, polyploidy, multiple NOR loci, and long generation times have all contributed to the concerted evolution of ITS sequences in mulberry.

Phylogenetic analyses of cpDNA sequences

Certain cpDNA sequences are widely used molecular markers in phylogenetic analyses (Nepal & Ferguson, 2012; Xu et al., 2012). With the development of universal primers for cpDNA and chloroplast genome sequencing, increasing numbers of studies have been conducted based on cpDNA (Huang et al., 2014; Wu et al., 2014). Another characteristic of cpDNA is that it is maternally inherited. To date, six chloroplast genomes of mulberry have been reported, and they can provide comprehensive information about the evolution of the whole chloroplast genome (Chen et al., 2016a; Hu et al., 2014; Kong & Yang, 2016, 2017; Ravi et al., 2007). Several other molecular markers have been used in phylogenetic studies, with trnL-trnF being the most commonly used (Ayinampudi et al., 2011; Nepal & Ferguson, 2012). In this study, more variable cpDNA regions (trnL-trnF and trnT-trnL) were used for phylogenetic analyses to study the concerted evolution of mulberry. The mulberry accessions in clade B were separated into two clades in the phylogenetic tree based on cpDNA sequences. The sequences in clade I were almost identical, sequences in clade II showed more variations and a subclade was clustered (Figs. 4 and 5). Successive evolution pattern of the sequences was shown in Fig. 4, suggesting a closer progenitor among them. Thus, clade II could be treated as a single clade. These results were consistent with those of previous studies in which two clades were proposed based on phylogenetic analyses of TrnL-TrnF (Zhao et al., 2005; Zhao et al., 2007).

Evolution of mulberry

Based on the analyses of ITS sequences, the 33 functional ITS sequences in clade B were separated into short-type sequence clades (clade 1 and clade 2) and a long-type sequence clade (clade 3). Two pseudogroups were also divided in the ITS phylogenetic tree. These findings suggest that there is a shallow level of phylogeny among 29 mulberry accessions, consistent with the two clades detected in the cpDNA phylogenetic tree (Bailey et al., 2003). Thus, clade B contains two species.

On the basis of a study of the genus Leucaena, it was proposed that some ITS pseudogene sequences are inherited from the maternal progenitor (Hughes, Bailey & Harris, 2002). In this study, we detected three connection patterns between ITS pseudogene sequences and cpDNA sequences (Figs. 1 and 4). Connection pattern (1) was between pseudogroup 2 and clade I, and was detected in Shanxitiansang, Xinyizhilai, Huosang, and Shimiansang. In this pattern, both the ITS pseudogene sequences and cpDNA sequences showed high sequence similarities. Connection pattern (2) was between pseudogroup 1 and clade II, and was detected in Yun7, Yun6muben, Baojing7, and Pisang2. The ITS pseudogene sequences of these four mulberry accessions displayed sequence variations, but the cpDNA sequences did not. Connection pattern (3) was between pseudogroup 2 and clade II, and was detected in Hanguodabaizhenzhu, Lunjiao109, and Agentingsang. In these three patterns, the ITS pseudogene sequences have connections with genetic information from the female progenitor. This provides more clues to trace the evolution process of ITS sequences and the Morus genus.

The three connection patterns and the incongruent distribution patterns between the cpDNA and the ITS phylogenetic trees imply that hybridization has contributed to the evolution of mulberry. Based on the discussion above, a simple evolutionary pattern is proposed for mulberry (Fig. 6). First, mulberry with short-type ITS sequences (^aAA) existed at a certain time and showed stronger ecological potential. Those resources then hybridized with mulberry with long-type sequences (^bBB). The ITS sequences gradually became short-type sequences through concerted evolution or hybridization. Finally, long-type ITS sequences (B) were left as pseudogene sequences (B) in the offspring, and most mulberry had short-type functional ITS sequences (^aAA, ^aAA^B, and ^bAA^B). Thus, mulberry with long-type ITS sequences (^bBB) are the original mulberry resources, and are being polluted by those with short-type ITS sequences during evolution. The actual evolutionary process of mulberry will be more complex than this evolutionary pattern, and dominated by hybridization.

In the phylogenetic trees constructed from ITS and cpDNA sequences (Figs. 3 and 5), all the clades contained accessions with different ploidy levels, and there were cross-links between the clades in different trees (e.g., Huai302 (2n = 12x = 84) belonged to clade 2 and clade I). These results indicate that hybridization and polyploidy have played important roles in the evolution of mulberry, even in the accessions with higher ploidy levels. Clade 2 of the phylogenetic tree constructed from ITS sequences contained mulberry accessions with different high ploidy levels, identical ITS sequences, and successively evolved cpDNA sequences, implying that concerted evolution has played a key role in the evolution of this clade. In conclusion, polyploidy, hybridization, and concerted evolution have all played important roles in the evolution of ITS sequences in mulberry.