cpSNP discovery and genotyping for a Pinus taeda breeding population with targeted comparison to related conifers

DNA extraction and sequencing

Fresh needles of P. taeda were collected from the Yingde Research Institute of Forestry in Guangdong Province, P.R. China (24°15’N, 113°25’E). Total cpDNA was extracted from approximately 100 g needles using a high-ionic-strength medium (Bookjans, Stummann & Henningsen, 1984). The integrity and purity of DNA were assessed by 1% agarose gel electrophoresis. The purity of the samples was further evaluated by measuring the OD260/280 ratio with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). The DNA concentration was quantified by an Invitrogen Qubit fluorometer. Libraries were prepared in accordance with the manufacturer’s instructions with an average insert size of 350 bp. Purified cpDNA was sequenced on an Illumina HiSeq 2500 platform (Science Corporation of Gene, Guangzhou, China) under a paired-end 100 bp mode. Raw read quality was then assessed using fastp (v0.23.2) with the following parameters: a cut window size of 8 nucleotides and a qualified quality Phred score threshold of 20.

De novo CP genome assembly

The chloroplast assembly software SPAdes (v3.15.5; Bankevich et al., 2012) was used to assemble the chloroplast genome of loblolly pine, with K-mer length set to 79 and 97. The assembly’s integrity was verified by mapping reads back to the consensus sequence using Burrows-Wheeler Aligner (BWA; v0.7.17-r1188) and SAMtools (v1.14). The P. taeda cp genome was compared with other available cp genomes of Pinaceae by using the CGView Comparison Tool (CCT; Grant & Stothard, 2008). Gene annotations were performed using the Clusters of Orthologous Groups (COGs) database. The Basic Local Alignment Search Tool (BLAST) was employed to align the other genomes to P. taeda, facilitating the identification of homologous sequences. The complete annotated genome was shown as a circular map using OrganellarGenomeDRAW (OGDRAW; Greiner, Lehwark & Bock, 2019). The adenine (A) and thymine (T) distributions were measured based on the results of the AT-skew equation, as follows: AT-skew = (A−T)/(A+T).

Codon usage

Codon usage, quantified as relative synonymous codon usage (RSCU) values (Sharp & Li, 1987), was determined for all protein-coding genes across the 33 Pinaceae species. Statistical analyses of RSCU distributions were performed, and results were visualized as a heatmap using the pheatmap package (v1.0.10) in R (R Core Team, 2022).

Ka/Ks analysis

The rates of nonsynonymous (Ka) and synonymous (Ks) substitutions within the chloroplast genomes of five Pinus species were estimated using seqinR (v3.4-5; Gouy et al., 1984; R Core Team, 2022). Histograms of these substitution rates were then created using ggplot2 (v3.0.0) to visualize the distribution of Ka and Ks values.

Single nucleotide polymorphisms analysis

There were 32 P. taeda cpDNA samples, selected from the Yingde Research Institute of Forestry in Guangdong Province, P.R. China, mixed and sequenced by the Illumina Hiseq 2500 platform with a paired-end 100 bp sequencing strategy. The resulting sequences were aligned against the known reference sequence of P. taeda (GenBank accession no. NC_021440.1). Subsequently, 18 single nucleotide polymorphisms (SNPs) were selected for the detection of 33 individuals in the breeding population of P. taeda. Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) was used for the primer design (Table 1), and the PCR products were sequenced using the Sanger method.

Phylogenetic analysis

To resolve the phylogenetic position of P. taeda, phylogenetic analyses were conducted using both complete cp genomes and the ycf1 gene sequences from 36 samples. Sequence data were retrieved from NCBI (https://www.ncbi.nlm.nih.gov/), comprising 24 Pinus species (P. rigida, P. jaliscana, P. oocarpa, P. elliottii, P. caribaea, P. taiwanensis, P. thunbergii, P. densata, P. wallichiana, P. strobiformis, P. bungeana, P. monophylla, P. nelsonii, plus P. taeda, P. greggii, P. contorta, P. sylvestris, P. tabuliformis, P. krempfii, P. strobus, P. massoniana, P. lambertiana, P. sibirica, and P. koraiensis represented by two samples each) and the outgroup Cedrus deodara.

Both datasets were processed identically: multiple sequence alignment used MAFFT v7 (Katoh, Rozewicki & Yamada, 2019), and Maximum Likelihood trees were inferred in MEGA X (Kumar et al., 2018) with 1,000 bootstrap replicates.

Topological congruence between whole-genome and ycf1-only trees was assessed by calculating the normalized Robinson-Foulds (RF) distance using ape (v5.8.1; Paradis & Schliep, 2019) in R v4.5.1 (R Core Team, 2022).

Genome comparison analysis

A total of 33,022,566 base pairs (bp) of raw reads with 38.41% GC content were generated with a paired-end 100 bp read length. The raw sequences of the P. taeda cp genome were deposited in GenBank with accession number PRJNA1159385, and the assembled sequences with accession number NC_021440.1. The size of the genome was 121,530 bp, with a GC content of 38.5% (Fig. 1). The positions of all the identified genes in the cp genome and their functional categorization are presented in Fig. 1.

The cp genomes of 14 species of Pinaceae (P. massoniana, P. lambertiana, Pseudotsuga sinensis var. wilsoniana, Cedrus deodara, Abies koreana, Larix sibirica, Pseudolarix amabilis, P. krempfii, Keteleeria davidiana, P. monophylla, Tsuga chinensis, Cathaya argyrophylla, Picea sitchensis, and Nothotsuga longibracteata) were selected for comparison with P. taeda by using CCT (Fig. 1). The sequence identity between P. taeda and other Pinaceae representatives was analyzed, revealing that certain regions were less conserved. For instance, ycf1 showed sequence identity lower than 70%, while ycf2 exhibited identities below 80%. Previous studies have confirmed that ycf1 has an important role in the evolution and classification of Pinus (Daniell et al., 2016; Handy et al., 2011; Georgolopoulos, Parducci & Drouzas, 2016).

Codon usage

The codon usage patterns among the 71 distinct cp protein-coding genes in P. taeda were examined (Fig. 2). The 71 protein-coding genes consisted of 20,255 codons, with the majority ending in A or T and the remainder ending in C or G. Codons terminating in A or T exhibited significantly higher RSCU values compared to those ending in C or G, as detailed in Table 1. Statistical analysis of codon usage distributions across 33 Pinaceae species revealed infrequent use of codons ending with G or C, as shown in the heatmaps (Fig. 2). Similar patterns have been reported in other cp genomes (Delannoy et al., 2011; Raubeson et al., 2007; Ying et al., 2016).

Among the 33 analyzed Pinaceae species, codon usage patterns were largely conserved, except for P. taiwanensis, P. thunbergii, and P. koraiensis. Phylogenetic clustering revealed that P. taiwanensis and P. thunbergii formed a distinct clade, divergent from other Pinus species. Additionally, Picea and Abies occupied phylogenetically distant positions both from Pinus and from each other. The clustering results were inconsistent with traditional classification results (Gernandt et al., 2005).

Ka/Ks analysis

The study analyzed the molecular evolutionary characteristics of seven chloroplast genes (clpP, ycf1, cemA, matK, atpA, atpB, and atpI) from five species belonging to Pinus. The Ka/Ks ratios of the 10 species pairs of these genes are shown in Fig. 3. Significant differences in the Ka/Ks ratios were found between the seven chloroplast genes. Ka/Ks >1 was detected in seven species pairs for clpP and in five species pairs for matK. For ycf1 and cemA, Ka/Ks >1 was detected in three species pairs each, while the remaining species pairs had Ka/Ks ratios close to 1. For atpA, atpB, and atpI, Ka/Ks >1 was detected in only one or two species pairs, with the remaining species pairs having Ka/Ks ratios close to 1.

SNP analysis of the chloroplast genome

Compared with the known reference sequence, 71 SNPs and four indels were identified from the 33 individuals in the breeding population of P. taeda (Table 2). Of these, 35 SNPs were located in the gene-coding regions of 11 genes, while the remaining SNPs were found in an intergenic region. The study found that 50 SNPs had mutation frequencies below 15%, which constituted 70.42% of the total number of SNPs. The lower mutation rate indicates that the cp genome of P. taeda was conservative. Notably, three SNP loci had mutation rates between 30% and 50%, with two located in gene-coding regions (ycf1 and ycf2). Specifically, ycf1 contained 17 mutations (23.94% of the total mutation sites), and ycf2 had three SNPs, including one with the highest mutation rate of 48%. R poB had three SNPs, two of which had the same mutation rate. R poC2 had two SNPs with mutation rates at 24% and 3%. R poC1 had two SNPs with a similar mutation rate, as did psaA. The intergenic regions between matK and chlB, psbJ and petA, and rps15 and psaC contained 14 SNPs, eight of which had mutation rates higher than 10%.

Individual distinction of the breeding population

A total of 18 SNPs were selected to detect 33 individuals from the breeding population of P. taeda. Of these, 13 SNPs were mapped to ycf1. Analysis revealed that 32 individuals exhibited SNPs across at 15 different loci, and 23 individuals could be distinguished based on these SNPs. Only one individual’s genome matched the reference genome, as shown in Table 3. Among the 23 distinguished individuals, one additional SNP was identified, and seven SNPs were found to be the most prevalent.

Furthermore, three additional SNPs were found in 14 individuals, accounting for 60.87% of the distinguishable SNPs. A total of 13 individuals had two SNPs, which represented 39.40% of the total number of individuals (33). Compared to previous studies, the cpDNA markers had a comparatively high level of variability (Walter & Epperson, 2010). The chloroplast simple sequence repeats (cpSSRs) of P. slvestris revealed a high level of intra-populational polymorphism across the entire range (Semerikov et al., 2014). Additionally, the variation among populations accounted for 99% of the total variance (Pazouki et al., 2016).

Additionally, the remaining 10 individuals that could not be distinguished, each of which had only one or two SNPs, could be classified into four groups. Introducing more SNPs could potentially help distinguish these individuals. Specifically, individuals P043, P054, and P064 exhibited identical SNP profiles. Given their common origin, they may descend from the same paternal lineage. Similarly, individuals G09, G10, and 243 also had the same SNPs. G09 and G10 had the same origin, while 243 was selected from a different stand. The loblolly pine in China, which were all introduced from the southern United States, might be related to each other.

Phylogenetic insights from cp Genomes and ycf1

Phylogenetic reconstruction using complete cp genomes resolved the 24 Pinus species and the outgroup C. deodara into two well-supported clades (Clade A and Clade B; Fig. 4A). The topology contained 33 nodes, with 26 nodes achieving maximal bootstrap support (100%), while the remaining nodes attaining values ranging from 64% to 99%. This bifurcation is consistent with the canonical classification proposed by Gernandt et al. (2005). Notably within Clade A, P. taeda and P. rigida converge into the same branch. P. taeda emerged as a focal lineage closely affiliated with P. jaliscana, P. oocarpa, P. greggii, P. elliottii, and P. caribaea, corroborating clade relationships reported by Duan et al. (2016) and Zeb et al. (2019); Zeb et al. (2022).

Analysis of ycf1 gene sequences from the same chloroplast genomic resources recovered an overall congruent phylogenetic architecture (Fig. 4B), despite minor topological variations. The resultant tree likewise recovered Clades A and B across 33 nodes, with 16 nodes exhibiting 100% bootstrap support and the remaining nodes having bootstrap values ranging from 52% to 99%. While P. jaliscana, P. oocarpa, and P. greggii retained their close clustering with P. taeda within Clade A, P. elliottii and P. caribaea failed to group with this lineage. Instead, they occupied distinct positions with bootstrap supports of 92% at their respective nodes.

Quantitative assessment of topological concordance revealed minimal divergence between the whole-plastome and ycf1-based phylogenies, evidenced by a normalized RF distance of 0.15. This strong congruence highlights the efficacy of ycf1 as a single genomic marker for resolving genus-level evolutionary relationships in Pinus.