Characterization of the complete chloroplast genomes of five Populus species from the western Sichuan plateau, southwest China: comparative and phylogenetic analyses

Plant materials and DNA extraction

The fresh leaves of P. cathayana were collected in Kangding (101°56′26″E, 29°59′36″N, Sichuan, China; altitude: 3,109 m), while the samples of P. kangdingensis, P. pseudoglauca, P. schneideri and P. xiangchengensis were collected in Kangding (101°36′43″E, 30°05′20″N, Sichuan, China; altitude: 3,554 m), Yajiang (100°54′06″E, 29°59′14″N, Sichuan, China; altitude: 3,598 m), Litang (101°36′43″E, 30°05′20″N, Sichuan, China; altitude: 4,018 m) and Xiangcheng (99°40′33″E, 28°55′47″N, Sichuan, China; altitude: 3,530 m), respectively. The voucher specimens of the five species were deposited at the herbarium of Southwest Forestry University, Kunming, China. Total genomic DNA was extracted with the Ezup plant genomic DNA prep kit (Sangon Biotech, Shanghai, China), and DNA samples were stored at −80 °C at the Key Laboratory of State Forestry Administration on Biodiversity Conservation in Southwest China, Southwest Forestry University, Kunming, China.

Genome sequencing, assembly and annotation

Total DNA was used to generate libraries with an average insert size of 400 bp, which were sequenced using the Illumina HiSeq X platform. Approximately 15.0 GB of raw data were generated from each genome with 150 bp paired-end read lengths. Then, the raw data were used to assemble the complete cp genome using GetOrganelle software (Jin et al., 2018) with P. trichocarpa as the reference. Genome annotation was performed with the program Geneious R8 (Biomatters Ltd, Auckland, New Zealand) by comparing the sequences with the cp genome of P. trichocarpa. The tRNA genes were further confirmed through online tRNAscan-SE web servers (Schattner, Brooks & Lowe, 2005). A gene map of the annotated Populus cp genome was drawn by OGdraw online (Lohse et al., 2013).

Indices of codon usage

As an important indicator of codon usage bias, the relative synonymous codon usage (RSCU) value is the frequency observed for a codon divided by its expected frequency (Sharp, Tuohy & Mosurski, 1986; Sharp & Li, 1987). The amino acid compositions and RSCU values of the five Populus cp genomes were calculated using the CodonW program (Peden, 1999). Because short CDSs generally result in large estimation errors for codon usage, CDSs shorter than 300 bp in length were excluded in codon usage calculations to avoid sampling bias (Rosenberg, Subramanian & Kumar, 2003). Finally, 58 CDSs from the five cp genomes were analyzed in this study.

Genome comparison

To investigate divergence in cp genomes, identity across the whole cp genomes was visualized using the mVISTA viewer in Shuffle-LAGAN mode for the five species, with the P. xiangchengensis genome as the reference. MAFFT version 7 software (Katoh et al., 2005) was used to align the five plastome sequences, followed by adjustment with BioEdit. To elucidate the level of sequence variation, we then performed sliding window analysis to assess the pairwise variability (Pi) over the plastomes in DnaSP version 5 software (Librado & Rozas, 2009). The window length was set to 600 bp, and the step size was set to 200 bp. The SNP variation was detected using the “find variation” function in Geneious R8.

Identification of simple sequence repeats and long sequence repeats

Simple sequence repeats in five Populus cp genomes were detected using Microsatellite identification tool (MISA) (Thiel, Michalek & Varshney, 2003) with the minimal repeat number set to 12, 6, 5, 5, 5, and 5 for mono-, di-, tri-, tetra-, penta-, and hexa nucleotide sequences, respectively. We used the online REPuter software to identify and locate forward (F), reverse (R), complemented (C), and palindromic (P) repeats. The following settings for repeat identification were used: (1) Hamming distance equal to 3; (2) minimal repeat size was set to 30 bp; (3) maximum computed repeats was set to 90 bp (Kurtz et al., 2001).

Gene selective pressure analysis of five Populus plastomes

To examine variation in the evolutionary rates of cp genes, we calculated the nonsynonymous substitution rates (Ka), synonymous rates (Ks), and their ratios (Ka/Ks) using model averaging in the Ka_Ks Calculator program according to the LWL85 method (Yang & Bielawski, 2000; Zhang et al., 2006).

Phylogenetic analysis

To explore the genetic relationships among the five species of the Populus genus, a total of 17 complete cp genomes of Populus and five plastomes of Salix were obtained from GenBank, and Itoa orientalis and Idesia polycarpa were used as the outgroups (Table S1). To examine the phylogenetic utility of different regions, phylogenetic analyses were performed based on the following data: (1) the complete cp DNA sequences, (2) the large single copy (LSC) region, (3) the small single copy region (SSC), (4) one inverted repeat region (IR), (5) the LSC+SSC region, (6) the LSC+SSC+IR region, and (7) a set of 85 common protein coding genes. All of the datasets were aligned using MAFFT under default settings. jModelTest 2.0 (Darriba et al., 2012) was used to determine the best-fitting model for each dataset based on the Akaike information criterion. A maximum likelihood method for phylogenetic analysis was performed based on the GTR+I+G model in RAxML version 8 (Stamatakis, 2014).

Features of the five Populus plastomes

The complete cp genomes of the five Populus species ranged from 156,465 (P. xiangchengensis) to 156,789 bp (P. cathayana) in length. The plastome size of P. schneideri was only one bp larger than that of P. pseudoglauca. The plastome size of P. kangdingensis was 11 bp larger than that of P. pseudoglauca and 166 bp smaller than that of P. cathayana (Table 1; Fig. 2). The five cp genomes included a pair of IRs of 27,620 bp in the three species P. kangdingensis, P. pseudoglauca and P. schneideri and an IR pair of 27,672 bp in P. cathayana and 27,570 bp in P. xiangchengensis. The GC contents were consistent in P. kangdingensis, P. pseudoglauca, and P. schneideri, with 34.5%, 30.5% and 42.0% in the LSC, small short SSC and IR regions, respectively (Tables 1 and 2).

Each of the P. cathayana, P. kangdingensis, P. pseudoglauca, P. schneideri, and P. xiangchengensis cp genomes encoded 130 functional genes; 112 of these were unique genes, including 78 protein-coding genes, 30 tRNA genes and four rRNA genes (Table S2). Most of these genes occurred as a single copy, while 18 genes were double copies: seven protein-coding genes, seven tRNA genes and four rRNA genes. The LSC region contained 59 protein-coding genes and 22 tRNA genes, whereas the SSC region contained 10 protein-coding genes and one tRNA gene.

Codon usage

Most protein-coding genes had the standard AUG sequence as the start codon, but ndhD started with GUG, and rpl16 started with ATC. ATC as an initiation codon has been reported in other cp genomes (Raubeson et al., 2007; Wu et al., 2017). GUG start codons have been reported in tobacco, but they are very rare in eukaryotic genomes (Kuroda et al., 2007). When GUG was the start codon of a protein, it was still translated as Met because of the separate tRNA used for initiation. Furthermore, the codon usage patterns of the 58 distinct protein-coding genes in the five plastomes were examined, and the plastomes of P. kangdingensis, P. pseudoglauca, and P. schneideri were consistent, with a length of 75,990 bp and encoding 25,330 codons, while those of P. cathayana and P. xiangchengensis were 75,864 and 75,840 bp in size and encoded 25,288 and 25,280 codons, respectively, as presented in Table S3. Coding ending with A and T/U had RSCU values >1 for the five Populus cp genomes, indicating that they were used more frequently than synonymous codons and may play major roles in the A+T bias of entire cp genomes. There was a general excess of A- and U-ending codons. All three stop codons were present, with UAA being the most frequently used among the five plastomes (Table S3). In addition, leucine (Leu, 10.67%, 10.65%, 10.65%, 10.65%, and 10.65%) and cysteine (Cys, 1.14%) were the most and least commonly coded amino acids, respectively, among the five plastomes (Table S3; Fig. 3).

Comparative analysis of the five Populus plastomes

In this study, the cp genomes of the five Populus species were well conserved, and no gene organization rearrangement occurred when P. xiangchengensis was used as a reference (Figs. 4 and 5). LSC, SSC, and IR sections of the three Populus species of P. kangdingensis, P. pseudoglauca, and P. schneideri were highly conserved and smaller than those of P. cathayana, while the IR regions were larger than those of P. xiangchengensis. Detailed comparisons of the IR-SSC and IR-LSC boundaries among the cp genomes of the five species are presented in Fig. 6. Two complete or fragmented copies of rpl22 and ycf1 were located at the boundaries between the LSC or SSC regions and IR regions among the five Populus plastomes. The rpl22 gene crossed the IR-LSC with only one bp variation in sequence length among the five plastomes. Parts of the ycf1 gene (15 (P. cathayana) -158 bp (P. kangdingensis, P. pseudoglauca, and P. schneideri)) were found in the SSC region at the IRb-SSC junction, whereas a portion of the ycf1 gene (1,689 (P. xiangchengensis) to 1,707 bp (P. cathayana)) was present in both IRs. A 61 bp overlap between ycf1 and ndhF was found in P. kangdingensis, P. pseudoglauca, and P. schneideri. The Pi values within the slide window of 600 bp in the five plastomes varied from 0.00001 to 0.00335 (Table 3), with a mean of 0.00210. However, nine highly variable loci (Pi > 0.01), including trnG-atpA, psbZ-trnfM, trnL-ndhJ, ndhC-trnV, ycf4-cemA, trnN-trnR, ycf1, ccsA+ccsA-ndhD, and trnR-trnN, were located in the five Populus plastomes (Fig. 7). Among these regions, trnG-atpA, psbZ-trnfM, trnL-ndhJ, ndhC-trnV, and ycf4-cemA were located in the LSC region, ycf1 and ccsA+ccsA-ndhD were in the SSC region, and trnN-trnR and trnR-trnN were in the IR regions.

Number and forms of mutations

We investigated SNPs, the most abundant type of mutation, in the five plastomes, with P. xiangchengensis as the reference. In gene-coding regions, we detected 70 SNPs in the comparative combination of P. cathayana–P. xiangchengensis, including 33 transition (Ts) and 37 transversion (Tv) SNPs, as well as 160 (97 Ts and 63 Tv), 166 (101 Ts and 65 Tv) and 164 (99 Ts and 65 Tv) SNPs in the plastomes of P. kangdingensis–P. xiangchengensis, P. pseudoglauca–P. xiangchengensis, and P. schneideri–P. xiangchengensis, respectively (Table 4). Furthermore, 106 (38 Ts and 68 Tv), 323 (130 Ts and 193 Tv), 316 (130 Ts and 186 Tv), and 314 (131 Ts and 183 Tv) SNPs were detected in noncoding regions among the four comparative combinations, respectively (Table S4).

In our study, a total of six small inversions (petA-psbJ, ndhC-trnV, trnN-trnR, ccsA-ndhD, ndhD-psaC, and ndhF-trnL) were identified based on the sequence alignment of the five complete cp genomes (Fig. 8). The small inversions from ndhC-trnV and ndhD-psaC occurred in only P. xiangchengensis, those from ndhF-trnL occurred in P. pseudoglauca and P. schneideri, those from trnN-trnR occurred in P. kangdingensis, P. pseudoglauca and P. schneideri, and those from ccsA-ndhD occurred in the four species other than P. cathayana, while the inversion from petA-psbJ occurred in the four species other than P. xiangchengensis (Fig. 8).

Synonymous (Ks) and nonsynonymous (Ka) substitution rate analysis

In this study, ratios of nonsynonymous (Ka) vs synonymous (Ks) substitutions were calculated for 78 shared unique protein coding genes in P. cathayana, P. kangdingensis, P. pseudoglauca, and P. schneideri, with P. xiangchengensis as the reference. Among these genes, only 19 protein-coding genes had Ka/Ks values (Fig. 9; Table S5). The Ka/Ks values of the remaining protein-coding genes could not be calculated because Ka or Ks was equal to 0, indicating that these sequences were conserved without nonsynonymous or synonymous nucleotide substitution. The Ka/Ks ratios of all genes except rpoC2 in P. pseudoglauca (1.00903) and the rbcL gene in P. kangdingensis (2.26407), P. pseudoglauca (2.26407), and P. schneideri (2.26407) were less than 1 (Fig. 9).

SSR and long repeat analysis

With MISA, a total of 170 SSR loci were detected, of which mononucleotide repeats (P1) comprised 148 (87.06%) of all SSRs and all of the mononucleotides composed of poly A (polycytosine) and poly T (polythymine) repeats (Table 5). Within the five plastomes, SSR loci were primarily located in the LSC region, followed by the SSC region. A total of 15 SSR loci were detected in the protein-coding genes rpoB, rpoC2, and rps8, with all others situated in intergenic spacers and introns (Table S6). A total of 28, 39, 39, 39, and 25 SSR loci were detected in the P. cathayana, P. kangdingensis, P. pseudoglauca, P. schneideri, and P. xiangchengensis cp genomes, respectively (Table 5). The corresponding numbers of these repeats in P. kangdingensis, P. pseudoglauca, and P. schneideri matched each other and consisted of 33 P1, two dinucleotide (P2) and four compound (C) repeats. Comparison among the five plastomes revealed that four P1 loci were found in only P. xiangchengensis, one C locus was found only in P. cathayana, and 15 SSR loci (11 P1, 1 P2, and 3C repeats) were detected in the plastomes of P. kangdingensis, P. pseudoglauca, and P. schneideri (Table S6).

In the plastomes of the five Populus species, we found 58 repeats in P. cathayana, which was a higher number than those found in the other four species (49, 48, 48, and 48 repeats, respectively) (Fig. 10A). P. pseudoglauca and P. schneideri shared the same number and types of repeats (17 forward repeats (F), 23 palindromic repeats (P), five reverse repeats (R), and three complement repeats (C)) (Fig. 10A). The majority of repeats (84.86%) varied from 30 to 39 bp in length (Fig. 10B). Variation in the number of repeat sequences situated in the four parts of the plastome was observed between species (Table S7).

Phylogenetic analysis based on the cp genome

Seven data partitions (complete cp genomes, LSC, SSC, IR, LSC+SSC, LSC+SSC+IR region, and protein coding regions) from 22 Salicaceae cp genomes were used to construct phylogenetic trees (Figs. S1–S6; Fig. 11). However, the best resolution in phylogenetic relationships was achieved using full cp genome sequences; thus, we discuss phylogenetic relationships mainly based on Fig. 11. All of Populus was divided into four main highly supported clades (Fig. 11). Three species of the section Turanga were clade I members. Clade II included seven species (P. adenopoda, P. alba, P. davidiana, P. qiongdaoensis, P. rotundifolia, P. tremula, and P. tremula × alba) in section Populus and one species in section Aigeiros (P. nigra). Clade III consisted of three species in section Tacamahaca (P. kangdingensis, P. schneideri, and P. yunnanensis) and two species in section Leucoides (P. lasiocarpa and P. pseudoglauca). Clade IV included the four species in section Tacamahaca (P. balsamifera, P. cathayana, P. trichocarpa, and P. xiangchengensis), one species in section Aigeiros (P. fremontii) and one species in section Leucoides (P. wilsonii). Our results showed that P. kangdingensis, P. pseudoglauca, and P. schneideri were in clade III, while P. xiangchengensis formed a sister relationship with 100% bootstrap support to P. cathayana in clade IV.