Distributed under Creative Commons Cc-by 4.0

The complete chloroplast genome sequence of Morus cathayana and Morus multicaulis, and comparative analysis within genus Morus L ABSTRACT Trees in the Morus genera belong to the Moraceae family. To better understand the species status of genus Morus and to provide information for studies on evolutionary biology within the genus, the complete chloroplast (cp) genomes of M. cathayana and M. multicaulis were sequenced. The plastomes of the two species are 159,265 bp and 159,103 bp, respectively, with corresponding 83 and 82 simple sequence repeats (SSRs). Similar to the SSRs of M. mongolica and M. indica cp genomes, more than 70% are mononucleotides, ten are in coding regions, and one exhibits nucleotide content polymorphism. Results for codon usage and relative synonymous codon usage show a strong bias towards NNA and NNT codons in the two cp genomes. Analysis of a plot of the effective number of codons (ENc) for five Morus spp. cp genomes showed that most genes follow the standard curve, but several genes have ENc values below the expected curve. The results indicate that both natural selection and mutational bias have contributed to the codon bias. Ten highly variable regions were identified among the five Morus spp. cp genomes, and 154 single-nucleotide polymorphism mutation events were accurately located in the gene coding region.


INTRODUCTION
Mulberry (genus Morus, family Moraceae) is widely distributed in Asia, Europe, North and South America, and Africa. The trees are cultivated in East, Central, and South Asia for domesticated silkworm and silk production, which is of economic importance. Some mulberry species are also valued for their hard wood, delicious fruit, bark for paper production, and multiple uses in traditional oriental medicine (Asano et al., 1994;Kim et al., 1999). Linnaeus classified the Moraceae family into order Urticales in 1753 using morphological characteristics. However, the Moraceae family has been repositioned into order Rosales on the basis of the phylogenetic relationships among some nuclear genes and chloroplast (cp) loci or genomes (Kong & Yang, 2016;Su et al., 2014;Zhang et al., 2011). Morus is a genus of trees in family Moraceae, and species classification within Morus is the subject of ongoing controversy. In the Annals of Mulberry Varieties in China (The Sericultural Research Institute Chinese Academy of Agricultural Sciences, 1993), more than 3,000 mulberry germplasm resources were classified into 15 species and four varieties. In the Flora of China (Zhou & Gilbert Michael, 2003), Morus is classified into two groups and 11 species.
With the development of sequencing technology in recent years, in addition to nuclear genome sequences, cp genes (matK and rbcl), gene spacer regions (trnH -psbA, atpB-rbcL, trnC-ycf6, and trnL-F), and cp genome information have been used to study plant molecular systematics (CBOL Plant Working Group, 2009;China Plant BOL Group et al., 2011;Kress & Erickson, 2007;Roy et al., 2010). A nuclear gene internal transcribed spacer (ITS) and three cp DNA fragments in a total of four sequences are generally recommended for investigation of near-edge species. In Morus, Zhao et al. (2005) used the ITS and trnL-F sequence to cluster 13 mulberry genotypes, representing nine species and three varieties, into five branches. There are no other reports of molecular systematics applicable to the differentiation and classification of Morus species.
Chloroplasts are photosynthetic organelles that occur widely in algae and plants. Although the cp genome is more conserved than the nuclear genome, many mutation events such as indels, substitutions, and inversions have been identified in cp DNA sequences (Ingvarsson, Ribstein & Taylor, 2003). As a DNA polymorphism index, indels and single-nucleotide polymorphisms (SNPs) can be used at a low taxonomic level. The DNA polymorphism rate between cp genomes was 0.15% for M. yunnanensis and M. balansae, 0.3% for M. indica and M. mongolica (Kong & Yang, 2016;Song et al., 2015), and 0.2% among four different Chinese ginseng strains (Zhao et al., 2015). These results indicate the presence of variable characteristics among the cp genomes of different species.
Complete cp sequences can provide insight into evolution and natural selection, and have been used in significant contributions concerning evolutionary mechanisms for species and chloroplasts (Asheesh & Vinay, 2012). To date, cp genomes have been reported for three Morus species: M. mongolica, M. indica, and M. notabilis (Chen et al., 2016;Kong & Yang, 2016;Ravi et al., 2006). However, there has been no research on codon biology or genome evolution. Here, we report the cp genomes for another two Morus species (M. cathayana and M. multicaulis) and compare the codon usage bias (CUB), sequence divergence, and mutation events among the five Morus spp. cp genomes (Table 1).

Analysis of simple sequence repeats (SSRs) and long repeat sequences
SSRs in M. cathayana and M. multicaulis cp genomes were detected using MISA (Thiel et al., 2003) with the minimal repeat number set to 10, 5, 4, 3, 3, and 3 for mono-, di-, tri-, tetra-, penta-, and hexa-nucleotides, respectively. The distribution, nucleotide composition, and polymorphism among SSRs were investigated. For long repeat sequences in five Morus species, the program REPuter was used to assess the number and location of all four type of forward match (F), reverse match (R), complement match (C) and palindromic match (P) (Kurtz et al., 2001). The identity of the repeat was limited to greater than 90% and the size was no less than 25 bp, respectively.

Indices of codon usage
The amino acid composition and relative synonymous codon usage (RSCU) values for M. cathayana and M. multicaulis cp genomes were calculated using Mega 5.05 (Tamura et al., 2011). Then, the GC content of the first, second, and third codon positions (GC1, GC2, and GC3, respectively) and the overall GC content of the coding regions (GCc) were calculated manually. GC3S is the GC content at the third synonymously variable coding position excluding Met, Trp, and three stop codons. The CodonW program on the Mobyle server (http://www.mobyle.fr) was used to calculate the GC3s and the effective number of codons (ENc). An ENc plot (ENc vs GC3s) was also analyzed. To identify differences in coding sequences and corresponding amino acids among the five Morus cp genomes, pairwise alignment of the nucleotide sequence in coding region of five sequences was performed using Clustal X 1.83 (Supplemental Information 1). SNP events and transition (Ts) or transversion (Tv) in the plastomes were counted and positioned. In addition, SNPs were further classified into synonymous (S) and non-synonymous (N) substitutions using Mega 5.05 (Tamura et al., 2011). Then, a Z test (P < 0.05) was carried out by bootstrap method (1,000 replicates). The ancestral states were inferred using the ML method. The gene classification was according to Chang et al. (2006).

Assembly and features of cp genomes
The  Table 1).
The M. cathayana and M. multicaulis cp genomes both encode 132 predicted functional genes, of which 112 are unique genes, including 78 protein-coding genes (PCGs), 30 tRNA and four rRNA genes with 63,873 bp, 2,208 and 4,524 bp, respectively; 18 were duplicated in the IR region, including seven PCGs and seven tRNA and all rRNA genes (Fig. 1). Fifteen genes (10 PCGs and five tRNA genes) contain one intron, and three PCGs (clpP, ycf3, and rps12) have two introns. As found in most other green terrestrial plants, the maturase K (matK ) gene in the M. cathayana and M. multicaulis cp genomes is located within the trnK intron; and rps12 is a trans-spliced gene; the 5 -end exon is in the LSC region and the other two reside in the IR region separated by an intron. The four rRNA genes and two tRNA genes of trnI and trnA are clustered as 16S-trnI -trnA-23S-4.5S-5S in the IR region, the same as in the cp genome of M. mongolica and M. indica (Kong & Yang, 2016;Ravi et al., 2006) and most of other plants (Mardanov et al., 2008;Wang, Shi & Gao, 2013;Wu et al., 2014).

Analyses of repetitive sequences
A total of 83 SSR loci, accounting for 973 bp, and 82 SSR loci, representing 949 bp in length, were detected in M. cathayana and M. multicaulis cp genomes (Table 2). Among these, there are 59 mono-, 7 di-, 3 tri-, 11 tetra-, and 3 penta-nucleotide repeats, respectively, in the M. cathayana cp genome; the corresponding numbers of these repeats in M. multicaulis are 63, 6, 2, 9, and 2. Mono-nucleotide repeats accounted for 71.1% and 76.8% of total SSRs in M. cathayana and M. multicaulis, respectively, similar to the levels found for M.    et al., 2015). Six SSRs, including one di-nucleotide of (TC) 6 , one tri-nucleotide of (TTC) 4 , two tetra-nucleotides of (AAAG) 3 and (TTCT) 3 , and two penta-nucleotides of (AAGGA) 3 and (TTTCT) 3 in M. multicaulis or (ATTTC) 3 in M. cathayana, contain at least one C or G nucleotide, and the SSRs have a high AT content (97.5%). Most of the repeats are located in noncoding regions-i.e., intergenic spacers and introns-except for 10, which were found in seven coding genes of ycf1 (3) (Table 2). A total of 14 long repeat sequences that are longer than 25 bp were identified in the five Morus spp. cp genomes, including one reverse, six forward and seven palindromic matches. Most of these repeats are located in intergenic spacers or introns, except for one 30 bp palindromic repeat, which is in the trnS genes (Table 3). Among all the long repeats, only four was detected in all the five cp genomes, which means that these repeats might create a diversity of cp genomes, and provide valuable information for phylogeny of genus Morus (Yang et al., 2014).

Codon usage patterns in M. cathayana and M. multicaulis cp genomes, and comparison with other three Morus spp.
As an important indicator of CUB, the RSCU value is the frequency observed for a codon divided by the expected frequency (Sharp & Li, 1986). RSCU is close to 1.0 when all synonymous codons are used equally without any bias; RSCU is greater or less than 1 when synonymous codons are more or less frequent than expected (Gupta, Bhattacharyya  & Ghosh, 2004). Codons ending with A and T have RSCU > 1 for M. cathayana and M. multicaulis cp genomes, indicating that they are used more frequently than synonymous codons, and may play major roles in the A+T bias of the entire cp genome. In terms of nucleotide composition, the most frequent codons are composed of T or A or their combination (i.e., ATT for Ile, AAA for Lys, AAT for Asn, and TTT for Phe); the least frequent codons have a high GC content (i.e., UGC for Cys, CGC and CGG for Arg, ACG for Thr, GCG for Ala, and CCG and CCC for Pro) (Table S1). Further analysis of the base composition at each codon position in the coding region revealed that the average GC content for GC1, GC2, and GC3 is lower than the AT content in the five Morus spp. cp genomes, and that the GC3 content is lower than the GC2 and GC1 content significantly ( Table 4). The regular occurrence of A/T in the third codon position supports the previous finding that mutation towards A+T is a strong driving force for the cp genome. The effective number of codons (ENc) is widely used as a measure of CUB, with range in 20-61, for genes with extreme bias using only one codon per amino acid or no bias using synonymous codons equally (Wright, 1990). GC3s is indicator of the level of nucleotide composition bias (Ahmad et al., 2013;Wright, 1990), and ENc plots (ENc vs GC3 S ) are a

Table 5 Transitions (Ts) and transversions (Tv) in the chloroplast genome of five Morus species.
Ts Tv useful indicator of the factors affecting codon usage, as well as the force of mutation or other factors. Predicted values lie on or just below the expected curve when a gene's codon usage is constrained only by the G+C mutation bias, and values are considerably below the curve when codon usage is subject to selection for codon optimization (Wright, 1990). Similar to the ENc plot for the Asteraceae family (Nie et al., 2013), the majority of genes among the five Morus spp. plastomes follow the standard curve, which indicates that the codon bias was mainly caused by nucleotide composition bias at the third position (GC3s) ( Table 4, Fig. 2). In addition, there are several genes lying below the expected curve; this suggests that in addition to mutational bias, codon usage variation among genes can also be influenced by selection force (Sablok et al., 2011;Wright, 1990).

Sequence divergence among five Morus spp. cp genomes
Comparative analysis of sequence divergence among the five Morus spp. cp genomes revealed nucleotide variability (Pi) values in the range 0-0.28533 with average of 0.00432, indicating that the differences among the five cp genomes is small. However, nine intergenic regions (trnT -trnL, petN-psbM, psbZ -trnG, trnG-trnfM, rps4-trnT, psbE-petL, rpl32-trnL,   trnL-trnF and ccsA-ndhD) and one intron (trnL) are highly variable, with much higher values (Pi > 0.008) than other regions (Fig. 3). All the 10 loci are in single-copy regions, not IR regions. These regions with highly variable loci are not random but are clustered in 'hot spots ' (Shaw et al., 2007;Worberg et al., 2007). It has been reported that rpl32-trnL is particularly highly variable in Machilus and Solanum plastomes, and in spermatophyte (Dong et al., 2012;Dong et al., 2015;Sarkinen & George, 2013;Song et al., 2015), and that trnL-trnF, trnT -trnL, rps4-trnT, and the trnL intron are highly variable between two species of genus Morus L. (Kong & Yang, 2016). The petN-psbM intergenic region was reported as part of the trnC-trnD intergenic region, and has been used in lower-level phylogenetic studies in flowering plants (Lee & Wen, 2004). The petN-psbM region and the clpP intron were successfully used in a study of phylogenetic relationships in peach species (Dong et al., 2012). In the present study, four rarely reported highly variable loci, psbZ -trnG, trnG-trnfM, psbE-petL, and ccsA-ndhD, were found in Morus plastomes. Phylogenetic studies at the species level in Morus have not been carried out using highly variable regions, and our analysis provides a basis for further phylogenetic study of Morus using these highly variable regions.

Numbers and pattern of SNP mutations
As the most abundant type of mutation, we investigated SNPs in gene coding regions of five Morus spp. cp genomes. There were 154 SNPs detected, including 95 Ts, 59 Tv, and 74 N and 80 S substitutions (Tables 5 and 6). Among the Tv, 45 are GC content changes, 14 are between A and T, and between G and C ( Table 5). The Kn, Ks and their ratio are indicators of the rates of evolution and natural selection (Yang & Nielsen, 2000). Relative to the reference ancestral states, 74 non-synonymous substitutions was observed in 28 of 79 protein-coding regions in five Morus spp. (Table 6). The photosynthetic apparatus gene psaB of M. indica share extra synonymous substitution sites, while ycf1, photosynthetic metabolism gene rbcL of M. notabilis and gene expression gene rps14 of M. indica has more non-synonymous than synonymous substitution sites, suggesting a relatively high evolution rate for these genes, as for M. yunnanensis and M. balansae (Song et al., 2015). A Z test can be used to detect positive selection by comparing the relative abundance of synonymous and non-synonymous substitutions within the gene sequences, the results of the Z test for each gene showed that rbcL and ycf1 genes of M. notabilis exist positive selection, and psaB of M. indica exists negative selection, which means that the genes of Morus chloroplast might exert different evolutionary pressures. For all of these substitutions, 74.03% and 25.32% of the SNP sites are in LSC and SSC regions, respectively, while only one SNP site is in M. indica IR region, consistent with the more conservative characteristics of IR compared to LSC and SSC regions (Dong et al., 2014;Jo et al., 2011).

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by the Youth Science and Technology New Star Program of Shaanxi Province, China (No. 2013KJXX-96). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.