The complete chloroplast genome sequences of three Broussonetia species and comparative analysis within the Moraceae

Sample collection, DNA extraction, and genome sequencing

Fresh leaves of B. papyrifera and B. kazinoki were sampled from a single plant of each species growing in Ankang (108°58′55′E, 32°41′50′N), Shaanxi Province, China, on September, 2015. The fresh leaves of B. kaempferi were sampled from a single tree growing in Liuzhou (109°26′59′E, 24°17′12′N), Guangxi Province, China, on April, 2017. The voucher specimens of the three species were planted in the Sericultural Research Base of Ankang University (Fig. 1).

A modified high-salt method was used to extract the genomic DNA from the leaves (Shi et al., 2012). The resulting genomic DNA was then fragmented and indexed by barcoding. The paired-end libraries, with an insert size of ∼350 bp, were then constructed in accordance with the standard Illumina protocol and the sequencing was carried out on the HiSeq 2000 platform with 125-bp paired-end reads (Illumina Inc., San Diego, CA, USA).

Genome assembly and annotation

The chloroplast reads were isolated from the raw reads using the bowtie2 software (Langmead & Salzberg, 2012) with a very sensitive local model and cp genome of Ficus racemosa L. (Mao & Bi, 2016) used as a reference. The resulting reads were assembled using SOAPdenovo2 with Kmer = 63 (Luo et al., 2012). Then, all of the contigs were mapped to the reference cp genome using BLAT (Kent, 2002) to identify their position and direction. Sequences with ambiguous alignment were trimmed manually and were regarded as gaps. The gaps were filled by the consensus sequence, which was generated with the model implemented in MAQ software (Li, Ruan & Durbin, 2008). The process was repeated until the reference genome was fully supported. SAMtools (Li et al., 2009) was used to parse the depth of assembly sequences.

The preliminary annotations of B. papyrifera, B. kazinoki, and B. kaempferi cp genomes were performed using the online automatic annotator, DOGMA (Wyman, Jansen & Boore, 2004). The protein coding genes (PCGs) and rRNA genes were then verified using BLASTN searches (e-value cutoff = 1e−10) against other Moraceae family cp genomes to ensure accurate annotations (Chen et al., 2015). The start and stop codons, or intron and exon junctions, of each annotated PCG were manually compared with the cp genomes of F. racemosa and Morus mongolica Schn. using the check_annotations.py module (Jin et al., 2020). The tRNA genes were confirmed using tRNAscan-SE 1.21, which specified that mitochondrial/cp DNA was the source (Schattner, Brooks & Lowe, 2005).

The analysis of RPL22 deficiency or transfer in genus Broussonetia

The primers, F22 (5′-GCAAACCAAAGAGAATGATGAC-3′), and R22 (5′-CGAGCGTCTACCATTATACCTAC-3′), were designed for the amplification of the RPS3–RPS19 region, including the inter-genic (IG) regions of RPS3–RPL22 and RPL22–RPL19, and the RPL22 gene. They were used to identify RPL22 gene deficiency in the cp genome. The PCR amplifications of the RPS3–RPS19 region were carried out with the genomic DNA of B. papyrifera, B. kaempferi, and B. kazinoki, as templates. The RPS3–RPS19 region of hybrid B. kazinoki × B. papyrifera, and M. mongolica were amplified at the same time. The PCR products were then cloned and sequenced. The resulting sequences were then aligned using ClustalX1.83 software (Thompson et al., 1997), and the corresponding results were adjusted manually according to the RPS3, RPL22, and RPS19 gene borders.

The primers, rpl22f (5′-ATAACCCCGTCCTCGAGCTT-3′) and rpl22r (5′-AGAAGAGAAGGACCAAGCGA-3′), were designed inside the RPL22 gene, according to F. racemosa and M. mongolica, and were used to identify whether the RPL22 gene was transferred to the nuclear genome. The RPL22 genes of M. mongolica were amplified simultaneously as a control. The amino acid sequence of RPL22 from F. racemosa and M. mongolica were subjected to TBLASTN analysis of the whole-genome sequence of B. papyrifera and the raw data (Gertz et al., 2006).

Phylogenetic analysis

To illustrate the phylogenetic relationships among members of the family Moraceae, 32 complete cp genomes were downloaded from GenBank (Table S1). Three species Rosa chinensis Jacq., Rosa minutifolia Engelm. and Rosa rugosa Thunb. from the family Rosaceae were used as the outgroups. The sequences of 66 PCGs present in all the 38 species (35 plus three outgroup species) were extracted using Python 3.6. Then, each of the 66 PCGs were aligned, accounting for frame shifts and stop codons, by MACSE with default settings (Ranwez et al., 2018). Phylogenetic analyses were performed using the concatenated nucleotide sequences by both the Maximum-Likelihood (ML) and Bayesian Inference (BI) methods. The generalized time reversible (GTR) with invariable sites (+I) discrete gamma (+G) model was selected by model test applying the Akaike information criterion (AIC) and Bayesian information criterion (BIC) using jModelTest (Darriba et al., 2012). ML phylogenies were inferred using RAxML 8.2.12 software, the bootstrap analysis was performed with 1,000 replications (Stamatakis, 2014). BI phylogenies were inferred using MrBayes 3.2.7 (Ronquist et al., 2012), with Markov chain Monte Carlo (MCMC) algorithm of 1,000,000 generations, sampled every 1,000 generations until convergence. The first 25% of the trees were discarded as burn-in, while the remaining trees were used to generate the consensus tree.

Divergence time estimation

The RAxML tree and software MCMCTREE of PAML (Phylogenetic Analysis of Maximum Likelihood) (version 4.9j) was used for divergence time analysis (Yang, 2007). We estimated the divergence time under the relaxed clock and the Hasegawa-Kishino-Yano 1985 (HKY85) nucleotide substitution model. The nucleotide substitution rate was set as (r = 1. 368 ×10⁻⁹, rgene gamma =2, 14.6, 1) gamma distribution (Muse, 2000; Xu et al., 2012). The two primary calibration points in our analyses were: (1) the divergence between the Moraceae and the Rosaceae, 42–161 million years ago (Mya), according to http://www.timetree.org; and (2) the divergence between Morus alba L. and Morus notabilis Schn., 6–17 Mya, based on genome sequence estimates (Jiao et al., 2020).

The MCMC method was run for 2,000 generations as burn-in, then sampled every ten generations until a total of 20,000 samples had been generated. Convergence at each node was determined using Tracer v1.7 by confirming the effective sample sizes (ESS) above 200, with the 95% highest probability density (HPD) accepted (Rambaut et al., 2018).

Natural selection event analysis

The natural selection events can be measured by non-synonymous substitutions (dN) and synonymous substitutions (dS). The observation dN >dS suggest positive selection, otherwise it will be negative selection. The value of dN/dS was calculated by two methods: (1) using HyPhy 2.2.4 software, the unrooted phylogenetic tree, the Branch-site model, and the Muse-Gaut 1994 (MG94) codon substitution model (Pond, Frost & Muse, 2005); and (2) using program CodeML in PAML and the pairwise comparison matrix of 35 concatenated nucleotide sequences (Yang, 2007).

PCG divergence and indices of codon usage

The analysis for the sequence divergence (pi) of each PCG among the 35 cp genomes was conducted using DnaSP v6 software (Rozas et al., 2017). The amino acid composition and relative synonymous codon usage (RSCU) values were calculated using Mega 11 (Tamura, Stecher & Kumar, 2021). The effective number of codons (ENc) is widely used as a measure of codon usage bias (CUB). GC3S indicated the GC content at the third synonymously-variable coding position and excluded Met, Trp, and the three stop codons, which are indicators of the level of nucleotide composition bias (Wright, 1990; Ahmad et al., 2013). ENc-plot (ENc vs GC3S) is a useful indicator of the factors affecting codon usage (Comeron & Aguadé, 1998; Peden, 1999). The values of ENc and GC3S were calculated using the CodonW program in the Mobyle server (https://mobyle.rpbs.univ-paris-diderot.fr/cgi-bin/portal.py#welcome). The values of GC12 (the GC contents at the first and second positions) and GC3 (GC contents at the third position) were also obtained by this method.

Genome assembly and features

The mean depth of the final assembled B. papyrifera, B. kazinoki, and B. kaempferi cp genomes were approximately 484.2-fold, 295.1-fold, and 608.2-fold, respectively (Table S2). Their lengths were 160,239 bp, 160,841 bp, and 160,592 bp, respectively (Table 1). The lengths were all shorter than those of Broussonetia kurzii Hook and Broussonetia luzonica Blanco, but longer than those of the genus Morus. All three cp genomes had the quadripartite structures typical of land plants, including a large single-copy (LSC) region, a small single-copy (SSC) region, and two inverted-repeat (IR) regions.

The GC contents of B. papyrifera, B. kazinoki, and B. kaempferi were 35.83%, 35.73%, and 35.64%, respectively. As with other members of the Moraceae, the GC content distribution of the three Broussonetia cp genomes was also uneven, being highest in the IR region, intermediate in the LSC region, and lowest in the SSC region (Table 1).

Overall, 131 genes were found in each of the three Broussonetia cp genomes, namely 86 PCGs, 36 tRNAs, and eight rRNAs (Table S3), arranged in the same linear order as in other species of the Moraceae, except for gene RPL22, which was absent from the three Broussonetia cp genomes (Fig. 2). Of all these genes, seven PCGs, seven tRNAs, and four rRNAs were duplicated in the IR regions. Additionally, ten PCGs (atpF, ndhA, ndhB, petB, petD, RPL2, RPL16, rpoC1, RPS16, and ycf68) and six tRNAs (trnI-GAU, trnG-UCC, trnK-UUU, trnL-UAA, trnV-UAC, and trnA-UGC) contained one intron, whereas three PCGs (RPS12, ycf3, and clpP) had two introns (Table S3).

The absence of the RPL22 gene from the genus Broussonetia

Agarose gel electrophoresis showed that the PCR amplification products of the F22 and R22 primers were different lengths (Fig. 3A). The sequences of the clones were the same as those in the HiSeq 2000 platform, which confirmed that the RPL22 gene is a truncated pseudogene in B. papyrifera and was lost in B.kaempferi, B.kazinoki, and the hybrid B. kazinoki × B. papyrifera chloroplast genome (Fig. 3B).

The PCR amplification products of the rpl22f and rpl22r primers contained no expected products, indicating that there was no transfer of RPL22 in B. papyrifera, B. kaempferi, B. kazinoki, and hybrid B. kazinoki × B. papyrifera (Fig. 3C). The results of TBLASTN analysis showed that no RPL22 gene was detected in the B. papyrifera genome, too. Also, further gene content analysis showed that RPL22 gene is lacked in other two Broussonetia plants (B. luzonica and B. kurzii), and Malaisia scandens Lour., a plant of genus Malaisia family Moraceae.

Phylogenetic and molecular dating analysis

The concatenated sequence, including 66 PCGs that were 56,217 bp in length and had 6,847 diverse loci, was used to construct the phylogenetic tree. Both ML and BI phylogenetic trees had two clades: clade 1 included Morus and Artocarpus, whereas clade 2 included the other seven genera. The three Broussonetia species in this study were clustered into one branch, and M. scandens was aggregated within this branch. B. luzonica and B. kurzii were present in the sister branch (Fig. 4). The genus Broussonetia was the earlier diverging lineage than the genera Ficus, Trophis and Antiaris, later than genus Streblus. Fifteen genes (infA, psbL, ycf3, clpP, ndhF, ycf1, RPS2, RPL16, RPL32, RPS3, RPL22, ndhA, cemA, ycf68, and ycf15) were lost to differing degrees from the 35 cp genomes (Fig. 4).

The likelihood of the best-scored ML tree was 137,376 and the bootstrap support values in all the nodes between the genera were greater than 90%, reaching 100% in five cases (Fig. 4). At the species level, the values at the partial nodes (six out of 35) were less than 90%, especially in the species of genus Morus. The values at the nodes of (Morus cathayana Hemsl. (M. alba, Morus indica L.)) were 30 and 67, respectively. The BI tree showed similar topologies with the ML tree. The three clades of M. cathayana, (M. alba, Morus indica) and M. mongolica in the BI tree could not be effectively distinguish (Fig. 4), which may be caused by the interspecific free hybridization within the genus Morus.

We further evaluated the molecular clock based on the phylogenetic tree, HKY85 nucleotide substitution, and two calibration points, 42–161 Mya for the divergence between the Moraceae and the Rosaceae and 6–17 Mya for the divergence between M. alba and M. notabilis. The mean and 95% highest posterior density (HPD) divergence times were mapped onto the phylogenetic tree. The common ancestor of Moraceae occurred before 52.74 Mya, whereas divergence between genus took place within the range of 33.01–50.10 Mya among the genera. The divergence of Broussonetia took place around 26.11 Mya. This was much earlier than the divergence of Morus, Artocarpus, and Ficus, which diverged approximately 5.50–19.98 Mya (Fig. 5).

Natural selection event analysis

The dN/dS ratios of all branches and species calculated using the unrooted tree were less than one except for M. cathayana, a wild mulberry species unique to China, which had a ratio of 1.906 (Fig. 5). The dN/dS ratios based on the pairwise comparison matrix of concatenated nucleotide sequences showed that the highest average value for any genus was 0.519 for Morus, compared with 0.229 for Broussonetia, 0.241 for Ficus, and 0.413 for Artocarpus (Tables S4 –S7). The ratio between M. cathayana and M. mongolica was 2.705 (Table S4), which was the highest ratio occurring in the genus Morus.

PCG divergence and indices of codon usage

The analysis of sequence divergence among the PCGs of the 35 cp genomes showed that the pi values were in the range of 0.0005 (psbF)-−0.0419 (RPS16) with a mean of 0.0174, indicating small differences among the PCGs. The pi values of PCGs in the SSC region were all greater than 0.01, with a mean of 0.0226, which was much higher than those in the IR regions, all of which were below 0.01 (Fig. S1).

A total of 26,288–26,306 codons were identified in the cp genome of three Broussonetia species; these were used in RSCU analysis. The results showed that the most commonly used codons consisted of mostly A and U (i.e., AUU for Ile, AAA for Lys, UUU for Phe, and AAU for Asn). The codons that were infrequently used consisted of more G and C (i.e., GGC for Gly, CGC and CGG for Arg, CUG and CUC for Leu, GCG for Ala, and CCG for Pro). In addition, there were 30 codons with RSCU >1; among them, 29 codons ended with A or U and one codon (UUG for Leu, 1.20–1.21) ended with G (Table S8).

The distribution of the average ENc values of the PCGs in this study was 25.6 (RPL36)–61.0 (psbF). The distribution of the average GC3S value was 12.10%–34.14%. The distribution of the majority of the genes in the ENc-plots fell on or near the trend line, except for psbF (ENc =61), which was above the curve. In addition, some genes (petN, psbI, psbJ, RPL33, and RPL36) fell far below the expected line, suggesting that the codon usage bias could be affected by many factors (Fig. 6).

Further analysis of the base composition showed that the distribution of the average GC content for GC12 was 33.55%–56.08%, and for GC3 was 14.32%–36.73%. The distribution range for each was relatively narrow. The neutrality plots of GC12 versus GC3 for the 35 species showed that the distributions of the genes deviated from the diagonal line (slopes = 1) and the slopes of the trend lines ranged from 0.036 to 0.171. The correlation coefficient ranged from 0.032 to 0.151 (Table 2), which indicated that there were no correlations between GC12 and GC3 for the genes.