DNA extraction and sequencing

The high-quality genomic DNA was isolated from approximately 100 mg samples using the E.Z.N.A.^® HP Plant DNA Mini kit (Omega Bio-Tek, USA). DNA quality and concentration were assessed in a Qubit 3.0 Fluorometer (Thermo Fisher Scientific Inc., USA), and DNA integrity was evaluated by a 1.0% (W/V) agarose gel.

Approximately 1μg of high-quality genomic DNA from each sample was used for the whole genome sequencing (WGS). The high-throughput sequencing library of C. poilanei was constructed and sequenced using the Illumina NovaSeq 6000 platform (NovoGene, Beijing, China) following the standard procedure of the manufacturer. The MGIEasy DNA Rapid Library Prep Kit (cat.# 1000006985; MGI-Tech., China) was used to construct the sequencing library of C. illinoinensis (cv. ‘Sioux’), C. cathayensis, C. dabieshanensis, C. hunanensis, C. tonkinensis, C. kweichowensis and C. sinensis, and the paired-end 100 bp reads were generated on the BGISEQ-500 platform (BGI, Shenzhen, China) according to the manufacturers’ procedures, respectively. The raw sequence data was submitted to the Sequence Read Archive (SRA) of NCBI ( Supplementary Table 1 ).

Data processes, assembly, validation and annotation

To obtain the high-quality plastome assemblies of 19 Carya species, 2 Gb raw reads were randomly selected for each species from NovaSeq 6000 data (C. poilanei), BGISEQ-500 PE100 data (seven species that were sequenced in this study) or from Illumina HiSeq X-TEN data retrieved from the NCBI database (eleven previously sequenced species), respectively. The low-quality reads from all samples, i.e., reads having over 50% bases with a quality value below 12, adaptor only and more than 10% N bases, were filtered out by SOAPnuke software (Chen et al., 2018), separately. The trimmed clean reads including nuclear and organelle genome data were used to assemble the plastomes for each species. The plastomes of 19 Carya species were assembled via NOVOPlasty software (version 3.7), with ‘Seed_RUBP_cp.fasta’ provided by the software as the seed input and the published Juglans regia plastome (accession no. NC_028617.1) as reference (Dierckxsens et al., 2016; Peng et al., 2017). The default K-mer size (K-mer = 39) was firstly selected and multiple K-mers were used for a synchronous assembling test during the assembling to obtain the ideal assemblies.

To verify the correction of the plastome assemblies, PCR amplification and Sanger sequencing were performed in the six EA species and the pecan variety ‘Sioux’, to further confirm the SC-IR boundaries of the assembled sequences, as well as several special genomic regions among EA species. The sequences of the primers are listed in Supplementary Table 2 .

The plastome assemblies were annotated using the online program GeSeq – Annotation of Organellar Genomes (Michael et al., 2017). Initial annotation, putative starts, stops, and intron positions were predicted against the annotation of Juglans regia and C. illinoinensis (GenBank accession numbers: NC_028617.1 and NC_041449.1) plastomes. The final annotations were determined by integrating the GeSeq prediction and manual correction. The circular plastome maps of Carya species were constructed using the online visualization program Organellar Genome DRAW version 1.3.1 (OGDRAW) (Stephan et al., 2019). The newly generated complete plastome sequences were deposited in GeneBank (Accession numbers were listed in Supplementary Table 3 ).

Comparative analyses of plastome structure features

The online program IRscope (Amiryousefi et al., 2018) was used to visualize the boundaries of LSC-IRb, IRb-SSC, SSC-IRa, and IRa-LSC of plastome among the 19 Carya species, The IR expansion and contraction among these Carya species were also compared afterward. DnaSP v. 6.12.03 (Julio et al., 2017) was employed to analyze the nucleotide diversities, sequence polymorphisms, and relative rates of plastome sequence divergence in the eighteen Carya species. In order to calculate the synonymous (Ks) and non-synonymous (Ka) substitution rates and the nucleotide variance (π and θ), we extracted the same individual functional protein-coding exons and aligned them separately using BioEdit v. 7.0.9.0 (Hall, 1999). To obtain selection patterns in protein-coding genes, we calculated the Ka and Ks values of each protein coding gene between two plastomes in 19 Carya species using DnaSP v. 6.12.03, and divided them to take the average to get the Ka/Ks ratio of each gene.

The online tool MIcroSAtellite (MISA, http://pgrc.ipk-gatersleben.de/misa/) was used to discover simple sequence repeats (SSRs or microsatellites) in the eighteen plastomes with the following parameters: ten for mono-nucleotide motifs, six for di-nucleotide, five for tri-nucleotide, and three for tetra-, penta- and hexa-nucleotides, respectively (Sebastian et al., 2017). The repeat numbers in the regions of LSC, SSC, IRa, and IRb were counted.

Codon usage bias was calculated using MEGA X v. 10.1.8 (Kumar et al., 2018) and EMBOSS v.6.3.1 (Peter et al., 2000; Itaya et al., 2013). In general, many genes less than 300 bp in length are believed have no real functions and the encoded proteins may not be accurate and do not necessarily contain domains. Too many such genes may interfere with the accuracy of the results. Therefore, the protein-coding genes with more than 300 nucleotides in length were extracted according to the annotation file of each plastome.

RNA-editing sites of protein-coding genes in Carya plastomes were predicted using the online program Predictive RNA Editor for Plant cp genes (PREP-Cp) suite with a minimal editing score of 0.7 (Mower, 2009) http://prep.unl.edu/).

Structural variations among Carya plastomes

Annotation-based circular maps of the 19 Carya plastomes can be categorized into three distinct structures ( Figure 1 ). The major structure (Structure I) is shared by 16 Carya plastomes, with the same number of genes (131 genes, representing 113 unique genes), gene order, and translation direction ( Figure 1 ). Eighteen unique genes were duplicated in the IR regions of these 16 plastomes ( Table 2 ). Structure II is the most different from the major structure (Structure I) with respect to the direction of transcription of rpl36, which uses the complement strand as a transcriptional template in C. sinensis ( Figure 1 ). A sequence loss of ~2.3 kb in length in the IRb region of C. poilanei and C. tonkinensis plastomes leads to distinguished genome structure (Structure III) and decreased gene content (a total of 129 genes). The loss of 2.3 kb sequences results in a contracted IR region and the loss of tRNA genes, trnN-GUU, and trnR-ACG in the IRb region ( Figure 1 ; Table 2 ).

IR junction variations among Carya plastomes

Despite relative conservation of IR/SC boundaries in plant plastomes, contraction and expansion of the IR-SC border regions are common in the process of plastid evolution, which is the major source of variation in angiosperm-plant plastome length (Raubeson et al., 2007; Wang et al., 2018). We examined the fluctuation of IR-SC borders together with the adjacent genes in the 19 Carya plastomes. Of the IR-SC junctions, IRa/LSC junction has the most conserved border within a gene spacer of rpl2 to trnH-GUG, with contraction sizes of 34 - 80 bp in the IRa region and 15 - 48 bp in the LSC region in the Carya plastomes ( Figure 2A ). Complete comparisons of these junctions revealed six patterns: Patterns I to VI ( Figure 2A ). Among them, Patterns I to IV are present in the 7 EA species and Patterns V and VI are present in the 12 NA species ( Figure 2 ).

Open in a separate window

Figure 2

IR junction variation and its diversity among the 19 newly assembled Carya plastomes. (A), Comparison of IRb-LSC, LSC-IRa, IRa-SSC, and SSC-IRb boundaries. (B), The modern distribution and the diversity of IR-SC boundaries.

In details, Pattern I includes two species, C. cathayensis and C. kweichowensis, in which the LSC/IRb border is located within the rps19 gene with the expansion of 1 bp and 4 bp from rps19 of LSC to the IRb region, and the IRb/SSC border is located within the ycf1 fragment with expansion 1 bp and 12 bp from ycf1 fragment of IRb region to the SSC region ( Figure 2A ). Pattern II consists of C. dabieshanensis and C. hunanensis, in which the LSC/IRb border is within the gene spacer of rps19-rpl2 with 1 bp contraction of rps19 in the LSC region and 34 or 76 bp contraction of rpl2 in the IRb region. The IRb/SSC border within the gene spacer of the ycf1-fragment and ndhF that contracted 5 bp and 14 bp of the ycf1 fragment in the IRb region and 350 and 32 bp in the SSC region, and their SSC/IRa border was similar to Pattern I ( Figure 2A ). Pattern III includes only C. sinensis, which similarly has borders of LSC/IRb and SSC/IRa with the same plastomes as in Pattern II and IRb/SSC border with the same plastomes as in Pattern I. Pattern IV contains C. poilanei and C. tonkinensis plastomes, with a distinguished IR-SC boundary pattern: IRb/SSC border within the ndhF gene with a 6-bp expansion from the SSC region to the IRb region, and a SSC/IRa border located in the gene spacer between trnR-ACG and rrn5S with contractions of 178 bp and 77 bp for both sides, respectively ( Figure 2A ). Meanwhile, the LSC/IRb border of the plastomes in Pattern IV is similar to those in Pattern II ( Figure 2A ). Of the 12 NA Carya species, 11 have the highly conserved SC-IR boundaries and were assigned into Pattern V, which has the LSC/IRb border within rps19 with 2 bp or 4 bp expansion to the IRb region, the IRb/SSC border on the last base of the ycf1-fragment in IRb and SSC/IRa border within the ycf1 gene with 1,093 bp expansion to IRa region ( Figure 2A ). Pattern VI includes only C. ovata, with the same borders of LSC/IR and SSC/IRa, except for the IRb/SSC border, which was within the ycf1-fragment with 5-bp expansion to the IRa region ( Figure 2A ). These results indicated more conserved SC-IR boundary patterns among NA Carya plastomes than those in EA that exhibit highly diverse IR-SC boundary patterns ( Figure 2B ).

Simple-sequence repeats among Carya plastomes

Simple-sequence repeats (SSRs), also known as microsatellites, are short tandem repeat DNA sequences that consist of repeating 1-6 nucleotide motifs widely distributed throughout the plastomes, which are important molecular markers for analysis of genetic diversity and relationships between species (Yang et al., 2011; Jiao et al., 2012). We detected 1,652 SSRs in the 19 Carya plastomes and the numbers of SSRs varied among the species, with a range from 76 in C. aquatica and C. texana to 106 in C. poilanei and C. tonkinensis ( Figure 3A ; Supplementary Table 5 ). Of these detected SSRs, mononucleotide repeats were the most abundant SSR motifs, which accounted for approximately 85.5% of the total SSRs, and over 99% of the mononucleotide repeats were composed of A/T repeats ( Figure 3A ; Supplementary Figure 4 ). While the tri- to hexanucleotide SSRs were occupied 7.75% of the total ( Figure 3A ; Supplementary Figure 4 ). In addition, the total number of SSRs in each Carya plastome also revealed significantly different patterns among species in EA (93 - 106) and NA (76 - 82) Carya species, in which the mononucleotide repeats (especially for A/T repeats) are the major contributor ( Figure 3A ; Supplementary Table 5 ). Meanwhile, the numbers of SSRs of the plastomes can also be classified into six patterns, corresponding to those of IR-SC boundary patterns ( Figures 2 , 3A ).

Open in a separate window

Figure 3

Comparison on SSR numbers, codon number and nucleotide diversities of Carya plastomes. (A), Statistics of SSR number for each Carya plastome. (B), Statistics of codon number for each Carya plastome. (C), The paired nucleotide diversities of Carya plastomes from species in Eastern Asia (EA) or North America (NA), or between two species from Eastern Asia and North America (EA-NA).

Codon usage analysis

Codon distribution analysis of 52 protein-coding genes with a length over 300 bp for each gene sequence in the Carya plastomes revealed a similar codon usage pattern among the plastomes of Carya species and J. regia (Supplementary Table 6). The protein-coding sequences contain a total of 21,190 to 21,210 codons in the Carya plastomes and 21,265 in J. regia, including stop codons. The total number of codons varied among EA Carya, whereas 11 of the NA Carya have the same number of total codons (21,210), with the only exception of NA species in C. ovalis plastome with 21,208 codons ( Figure 3B ). The total number of codons for each species can also be classified into six patterns as mentioned above.

Other than the codon distribution, we also counted the average effective numbers of codon (ENC) in the plastomes. The results showed that 17 of the 19 assembled Carya plastomes had similar codon usage bias (ENC ranged from 48.883 to 49.001), which were slightly lower than that of J. regia (48.776) (Supplementary Table 6). C. tonkinensis and C. poilanei have ENC = 56.51 and 56.52, which are significantly higher than any of other species. All the Carya plastomes have similar GC content in the codons (37.32% ~ 37.39%) except for C. tonkinensis and C. poilanei plastomes which have lower GC contents for the first (39.59%) and second codons (31.12% and 31.18%) and relatively high GC content (41.26% and 41.20%) for the third position ( Supplementary Table 6 ). Among all codons in the 19 Carya plastomes, leucine is the most abundant amino acid, while cysteine showed the least abundance ( Supplementary Figure 6A ). The most frequent synonymous codon in the Carya plastomes was AUU, which encodes isoleucine, and the least frequent codon, except for the stop codons, is UGC, which encodes cysteine ( Supplementary Figures 5, 6B ). Of the 20 amino acids and stop codons, only leucine each has six codon types, while methionine and tryptophan preferred one codon type in the plastomes of 19 Carya species and J. regia ( Supplementary Figures 6 ). In general, AUG encodes not only methionine but also the universal start codon for the nuclear genome of eukaryotes (Thach et al., 1966). We detected that the genes of ndhD, rpl16, ycf15, and rps19 use ACG, ATC, CTG, and GTG encodes (not ATG) to encode the start codon – AUG respectively.

We also assessed the relative synonymous codon usage (RSCU), a good indicator for measuring nonuniform synonymous codon usage bias in coding sequences (Lee et al., 2010). Both methionine and tryptophan have RSCU = 1 (Supplementary Figure 7), indicating that codons AUG and UGG have no bias or preferences. The plastomes of 19 Carya species and J. regia had 30 biased codons with RSCU > 1, and their third position is A/U except for leucine (UUG). UUA has the highest RSCU values (1.94 to 2.03) in leucine in all the plastomes of Carya species and J. regia, and AGC has the lowest RSCU values in leucine (0.38 for Carya species and J. regia). Moreover, leucine showed A or T (U) bias in all synonymous codons: UUA, UUG, CUU, CUC, CUA, and CUG. We observed that the RSCU value for the specific amino acid increased with the number of codons. The following codons have high RSCU frequency (>30%) and fraction: GAU (aspartate), GAA (glutamate), AUU (isoleucine), AAA (lysine), UUA (leucine), AAU (asparagine), UAU (tyrosine), UUU (phenylalanine), and CAA (glutamine); and the bias toward these nine codons was consistent with the low content of GC in the third codon position ( Supplementary Figure 7 ). Our analysis also showed that the RSCU value increased with the quantity of codons coding for a specific amino acid.

Prediction of RNA editing

We found that 62 post-transcriptional RNAs have editing modifications in 23 protein-coding genes ( Supplementary Figure 9 ; Supplementary Table 7 ). Genes of ndhB (11 editing sites) and ndhD (10 editing sites) contain the most RNA editing sites, followed by rpoB and rpoC2, with 6 and 5 editing sites respectively. Only one to three RNA editing sites were found in the rest of the 19 protein-coding genes. All of RNA editing sites potentially caused the conversion from cytosine to uracil after transcription, and 43 of them (69.4%) took place at the second nucleotide of codons, with two times the conversion rate compared to the first nucleotide (19, 30.6%) ( Supplementary Figure 8 ). However, no correlation was observed between gene length and gene number ( Supplementary Figure 8 ). Approximately 77% of the RNA modifications resulted in the conversion of hydrophilic to hydrophobic amino acids, mainly serine to leucine or phenylalanine (S to L, 22 editing sites; S to F, 8 editing sites), or proline to leucine (P to L, 8 editing sites; Supplementary Table 7). We also observed the conversion from proline to serine in three editing sites, representing the changes of amino acids from nonpolar to polar. On the species level, the majority of the RNA-editing sites were shared by all the Carya plastomes and only a few were species-specific, for example, CUC to UUU (L to F) was found only in C. poilanei at the amino acid position 71, CGG to UGG (R to W) only in C. sinensis, CCA to UCA (P to S) only in C. floridana, and CUU to UUU (L to F) only in C. laciniosa (Supplementary Table 7).

Analyses of selection and nucleotide diversity

In order to obtain selection patterns in protein-coding genes, the nonsynonymous (Ka) and synonymous (Ks) substitutions ratios (Ka/Ks = ω) were determined for 79 unique protein-coding genes with the comparison of the 19 Carya plastomes ( Supplementary Figure 9 ; Supplementary Table 8 ). Among the protein-coding genes, ndhA and petA have Ka/Ks ratios greater than 1.0 ( Supplementary Figure 10 ; Supplementary Table 8 ), indicating the genes were subjected to positive selection. The rest of the 77 protein-coding genes have Ka/Ks < 1.0 ( Supplementary Figure 9 ; Supplementary Table 8 ), indicating that purification selection happened only during plastome evolution. Specially, seven genes (atpF, ccsA, matK, rbcL, rpoA, rps15 and ycf1) have the Ka/Ks ratios between 0.5 and 1.0 ( Supplementary Figure 9 ; Supplementary Table 8 ). These results clearly indicate that protein-coding genes in plastomes of different plant species were subjected to diverse selection pressures.

Two parameters, Pi (π) and theta (θ), were used for measuring the nucleotide variability among protein-coding genes of the 19 Carya plastomes. The value πvaried among the protein-coding genes with a range from 0 to 0.08495 (average π = 0.00326) (Supplementary Table 8). The gene rpl36 (π = 0.08495) was notably variable among the protein-coding genes in the 19 Carya plastomes. The gene ycf1 showed a relatively lowerπ value (0.00493) among the protein-coding genes in the Carya plastomes, although it was commonly used as a representative plant DNA barcoding region (Dong et al., 2015). The variation of the other parameter θ showed the same patterns as the π value among Carya plastomes of the 79 protein-coding genes ( Supplementary Table 8 ). In comparison, plastomes of EA Carya species revealed higher nucleotide diversities (π and θ) than those in NA Carya species ( Figure 3C ).

Complete plastomes: An excellent tool for inferring the matriarchal phylogeny and geography of Carya

Our maternal phylogenetic analyses showed that six of the BI-based and all ML-based topologies inferred by different datasets of plastomes support the sister relationship of Carya and the monophyletic Cyclocarya-Juglans-Platycary-Pterocarya clade in Juglandoideae ( Figure 4 ; Supplementary Figures 10–13 ). Our results also strongly support two monophyletic subclades corresponding to the disjunctive geographical distributions of Carya in EA and NA ( Figure 4 ; Supplementary Figures 10–13 ). The results are highly consistent with those of previous inferences on the intergeneric phylogenetic relationships in Juglandoideae and the relationships between two subclades in the genus Carya, based on the molecular markers from partial or complete nuclear or/and organelle sequences (Manos and Stone, 2001; Manos et al., 2007; Zhang et al., 2013; Zhang et al., 2019; Mu et al., 2020; Zhang et al., 2021; Zhou et al., 2021). It is known that the plastome is a key organelle in the plant cells, performing photosynthesis and other metabolic processes related to the adaptation of the plant to its environments (Dierckxsens et al., 2016). Although the basic structure of the plastome is highly conserved throughout land plant lineages, it has been proven that differences in the sizes of the complete genome and the protein-coding gene content of the different genome regions related to the difference in selection pressure are informative in phylogeny and evolution for many plant lineages. Although the matriarchal phylogenetic relationships among species and/or varieties within the EA and NA subclades varied among the topologies inferred from different datasets of plastomes, the plastome-based physiologies still provide very important clues to the backbone relationships between EA and NA Carya ( Figure 4 ; Supplementary Figures 10–13 ). The plastome-based phylogeny of Carya could also provide a good example for exploring the evolution of plastomes as we present in this study.

By comparing the PP and BI values of all nodes among six topologies, almost all the nodes of the BI tree reconstructed from full-length plastomes showed the highest support (PP > 0.8) ( Figure 4 ; Supplementary Figures 10–13 ). Therefore, the full-length plastome dataset-based BI tree is considered the main reference for further analyses of maternal phylogenetic relationships among Carya species. The maternal phylogenetic relationships among EA Carya species are highly in agreement with their geographical distributions ( Figures 2B and 4 ; Supplementary Figure 17 ). The maternal phylogenetic positions of the species in NA subclade provide strong support for their overlapping and adjacent distribution patterns ( Figures 2B and 4 ; Supplementary Figure 17 ), although the topology does not match the morphological features-based division of two sections (Manos and Stone, 2001; Huang et al., 2019). Meanwhile, the topologies among species within EA and NA subclades exhibited significant deviation from the pattern presented by nuclear genome data (Huang et al., 2019). Our results demonstrated that the different contributions of nuclear and organelle genomes, compared to the control, based on the morphological features associated with the different evolutionary rates between organelle and nuclear genomes, showed the difference between topologies based on the different datasets of plastomes in the present study. In general, the nonparental-inherited plastomes are highly conservative in protein-coding gene functions and genome structure, which can also reveal the maternal origin and diversity of a plant (Raubeson and Jansen, 2005). However, the morphological characteristics could be controlled by parental-inherited materials in the nucleus and may be influenced by environment and selection resulting in adaptation, convergent evolution, and species diversity (Givnish, 2016). Taking this into account, it is easy to understand the contradiction between topologies built from plastomes and nuclear data, as both had different evolutionary trajectories.

Carya poilanei, has been regarded as an extinct species in the subfamily Juglandoideae for more than 60 years before its rediscovery in the Ailao Mountain area, Jianshui County, Yunnan Province, China (Zhang et al., 2022). This species has historically been classified as a member of Juglans (Chevalier, 1941), then moved to the genus Carya (Leroy, 1950), and finally botanically characterized (Leroy, 1955). Although it is classified in the section Sinocarya of the genus Carya based on its morphological traits, the phylogenetic and maternal relationships are still unclear. The rediscovery of the species brings this issue back into consideration. Our results indicate that C. poilanei has the highest similarity with C. tonkinensis in plastome structure and features, and it has the closest maternal phylogenetic relationship with C. tonkinensis, which is supported by their closest geographical distribution ( Figures 1 – 4 ; Tables 1 , 2 ; Supplementary 17 ). Therefore, we assume that these two species historically shared a common maternal ancestor. In addition, the taxonomic placement of four contentious species C. dabieshanensis, C. (Annamocarya) sinensis, C. glabra, and C. ovalis has been disputed for decades (Chang and Lu, 1979; Thompson and Grauke, 1991; Manos and Stone, 2001; Grauke and Mendoza-Herrera, 2012; Grauke et al., 2016). C. dabieshanensis has been treated as a member of C. cathayensis because of their high similarities in morphological features, except for the larger fruit size of C. dabieshanensis (Liu and Li, 1984). Our previous nuclear genome analyses also supported that C. dabieshanensis is close to C. cathayensis, but not to any members in the section Sinocarya (Huang et al., 2019). This study provided a full comparison of the morphological features between these two species and found no additional differences beyond the genomes ( Supplementary Figure 17 ). However, the plastome features of both species reveal a significant variation in the patterns of IR-SC junction: C. cathayensis belongs to Pattern I with C. kweichowensis, but C. dabieshanensis belongs to Pattern II with C. hunanensis ( Figure 2A ). The phylogenetic topology built from the complete plastomes also supports that C. dabieshanensis has the close matrilineal relationship to C. hunanensis but not to C. cathayensis ( Figure 4 ). Although taxonomists distinguished C. (Annamocarya) sinensis from the genus Carya by its distinctive taxonomic, botanical, and horticultural characteristics, high throughout genome-wide sequencing technology has provided a credible molecular phylogeny to distinguish plant species (Favre et al., 2020). The plastomes-based phylogenetic tree here also provided solid support that C. (Annamocarya) sinensis is a member of EA Carya ( Figure 4 ). From the plastome phylogenetic tree, C. glabra and C. ovalis were separated in the NA Carya subclade and this is also confirmed by their plastome features, the patterns of IR-SC junction variations ( Figures 2A , 4 ), as well as their morphological differences (Grauke, 2003). As such, the matriarchal-originating plastomes provide an effective tool for taxonomic placement of the outcrossing plant species, at least at genus level.