The complete chloroplast genome sequence of Morus cathayana and Morus multicaulis, and comparative analysis within genus Morus L

Sample collection, DNA extraction, and sequencing

Samples of M. cathayana were collected from the Qinba Mountain Area and M. multicaulis was acquired from Shandong Sericultural Research Institute. Both plants used in the study are now maintained in the mulberry field of Shaanxi Key Laboratory of Sericulture, Ankang University, China. The cp DNA was extracted from 10 g of fresh leaves using a modified high-salt method (Shi et al., 2012) and treated according to a standard procedure (Bartram et al., 2011) for sequencing on an Illumina Hiseq 2000 platform (2 * 125 bp).

Chloroplast genome assembly and annotation

Both reads were assembled using SOAP de novo software (Luo et al., 2012) with the Kmer and maximal read length set to 60 and 50, respectively. The resulted contigs were aligned to the M. mongolica cp genome as the reference genome. Ambiguous sequences were manually trimmed. The mean fold coverage of the final assembled M. cathayana and M. multicaulis cpDNA reached approximately 690- and 770-fold, respectively. Transfer RNAs (tRNAs), ribosomal RNAs (rRNAs), and protein-coding genes (PCGs) in the two cp genomes were annotated using Dual Organellar GenoMe Annotator (DOGMA) software (Wyman, Jansen & Boore, 2004). The tRNAs and their corresponding structures were further verified and predicted using the tRNAscan-SE 1.21 program (Schattner, Brooks & Lowe, 2005). The physical maps were drawn using the web tool Organellar Genome DRAW (OGDraw) v1.2 (Lohse, Drechsel & Bock, 2007).

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.

Sequence divergence and mutation events analysis

A sliding window analysis was conducted using DnaSP 5.0 software (Librado & Rozas, 2009) for comparative analysis of the sequence divergence (Pi) among the cp genomes of M. cathayana, M. multicaulis, M. mongolica, M. indica, and M. notabilis. The window length was 600 bp and the step size was 200 bp.

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 complete cp genomes of M. cathayana and M. multicaulis are both closed circular molecules of 159,265 and 159,103 bp (GenBank accession number: KU981118, KU981119), respectively (Fig. 1). Both cp genomes show the typical quadripartite structure of most angiosperms, and comprise a pair of IRs (25,639 bp for M. cathayana and 25,677 bp for M. multicaulis) separated by the LSC (88,143 bp for M. cathayana and 87,940 bp for M. multicaulis) and SSC (19,844 bp for M. cathayana and 19,809 bp for M. multicaulis) regions. The GC content of M. cathayana and M. multicaulis cp genomes is 36.16% and 36.19%, respectively. Similar to values reported for Rosaceae cp genomes (Su et al., 2014; Wang, Shi & Gao, 2013), the GC content of the five Morus spp. cp genomes is 36.16–36.37%, with uneven distribution within the cp genome: the GC content is highest in the IR (42.89–42.95%), intermediate in the LSC (33.77–34.12%), and lowest in the SSC region (29.20–29.35%) (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).

Table 2:

SSRs in M. cathayana and M. multicaulis chloroplast genomes.

SSR	Size	NUMBER	LOCI
a	10	8	3950, 5052,590,29073, 50058, 68954, 68987, 114495(ndhF)
	11	4	, 54282${}^{\underset{¯}{\mathrm{a12}}}$, 63153, 87796, 116614${}^{\underset{¯}{\mathrm{a10}}}$
	12	2	13584${}^{\underset{¯}{\mathrm{a13}}}$, 84902
	13	1	128314${}^{\underset{¯}{\mathrm{a14}}}$
	14	1	74502
	15	1	9565${}^{\underset{¯}{\mathrm{a11}}}$
	16	2	4799${}^{\underset{¯}{\mathrm{a13}}}$8984${}^{\underset{¯}{\mathrm{a15}}}$
t	10	17	5231, 9782, ,24363, 30678, 30944, 54323, 55220, 57416(atpB), 62926, 67243, 69081, 71234, 74300, , 117117,122510, 30637(ycf1), 132394(ycf1)
	11	8	479, 8575${}^{\underset{¯}{\mathrm{t10}}}$,34221, 57867${}^{\underset{¯}{\mathrm{t12}}}$,59882, 75017${}^{\underset{¯}{\mathrm{t13}}}$,79022,131496(ycf1)
	12	8	12693, 13275${}^{\underset{¯}{\mathrm{t13}}}$, 14182${}^{\underset{¯}{\mathrm{t10}}}$,27623(rpoB), 68830${}^{\underset{¯}{\mathrm{t13}}}$, 69893${}^{\underset{¯}{\mathrm{t11}}}$,72813,86135
	13	3	9207${}^{\underset{¯}{\mathrm{t14}}}$, 52130${}^{\underset{¯}{\mathrm{t12}}}$,128735
	14	1	64181${}^{\underset{¯}{\mathrm{t13}}}$
	16	1	81690${}^{\underset{¯}{\mathrm{t13}}}$
	17	1	49756
	18(20)	1	87252${}^{\underset{¯}{\mathrm{t14}}}$,
at	5	2	69153,116007(ndhF)
	6	1	10796
ta	6	3	5495, 21235(rpoC2), 119074
tc	6	1	64908(cemA)
aat	4	1	128715
ttc	4	1	71261
tat	4	1	50114
aaag	3	1	135481
aaat	3	3	24057, 38611, 47006
atta	3	2	33937, 116796
attt	3	2	14173, 62456
tatt	3	1	24394
tctt	3	1	111916
ttat	3	1	118221
aagga	3	1	14008(atpF)
atata	3	1	117818
atttc	3	1	24297[tttct3]

DOI: 10.7717/peerj.3037/table-2

Notes:

Parentheses, containing coding regions; boxed type, absent in M. cathayana but present in M. multicaulis; bolded type, absent in M. multicaulis but present in M. cathayana; underline and superscript, nucleotide length polymorphism; bracket, nucleotide content polymorphism, others are identical.

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. mongolica and M. indica cp genomes, and the study on lower and higher plants (George et al., 2015). Six SSRs, including one di-nucleotide of (TC)₆, one tri-nucleotide of (TTC)₄, two tetra-nucleotides of (AAAG)₃ and (TTCT)₃, and two penta-nucleotides of (AAGGA)₃ and (TTTCT)₃ in M. multicaulis or (ATTTC)₃ 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), ndhF(2), rpoC2, rpoB, atpB, atpF and cemA. The 10 SSRs in coding genes also occur in the M. indica and M. mongolica cp genomes. Comparison between M. cathayana and M. multicaulis revealed that 57 loci are identical, 19 exhibit length polymorphisms, six SSRs loci only exist in M. cathayana, and five SSRs loci only exist in M. multicaulis. One nucleotide content polymorphic of (TTTCT)₃ and (ATTTC)₃ exhibits in M. multicaulis and M. cathayana, respectively (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 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).