Characterization of the chloroplast genome of a relict tree, Pterocarya fraxinifolia (Juglandaceae), and its comparative analysis

Relict species have always excited evolutionary biologists and biogeographers who consider these species ‘living fossils’ or relics of prehistoric periods^1,2. These species have great value as research models for the geographical distribution of intercontinental rifts and as species that ensure biodiversity and ecosystem balance. Relict species also provide relevant information about the adaptation of species to specific environmental changes, as well as the impact of climate change on the animal and plant kingdoms³.

Hyrcanian forests are hotspots for biodiversity and are home to numerous relict species⁴, including 280 endemic and subendemic species^5,6,7. The genus Pterocarya Kunth (Juglandaceae), commonly referred to as wingnuts, has a disjunct distribution in East Asia and the Caucasus region with its most recent common ancestor present 40 Ma⁸. Pterocarya comprises six species, which are classified into two sections, Pterocarya (P. fraxinifolia, P. hupehensis, P. stenoptera, and P. tonkinensis) and Platyptera (P. macroptera and P. rhoifolia), on the basis of the presence or absence of scales on the terminal buds⁹. P. fraxinifolia is the only species in western Asia¹⁰. The remaining species of Pterocarya occur in eastern Asia, such as China and Japan^{11,12,13,14,15}. Recently, a series of studies have focused on the phylogeny, biogeography, population genetics, and landscape genetics of species in this genus^14,15. However, resources regarding the chloroplast genome in this genus are insufficient, and more research is still needed.

Pterocarya fraxinifolia is a deciduous tree that can reach 20–25 m in height and 1.8 m in trunk diameter and is wind-pollinated to produce wing-nut fruits¹². This species is among the last surviving trees growing in small scattered populations in the southern Caucasus region, which includes northern Iran, Georgia, Armenia, Azerbaijan, and the Anatolian region in Turkey^12,15. However, less than two decades ago, small populations were first recorded in western Iran in the provinces of Lorestan and Ilam in the Zagros Mountains¹⁶.

The chloroplast genome is widely used in phylogenetic studies because of its relatively conserved structure^17,18 and uniparental inheritance^19,20,. Chloroplast genomes can provide important information about the adaptation of species to different environmental conditions^21,22,23. Despite the slow evolutionary rates of chloroplast genomes, coding and noncoding regions are useful for the identification of closely related species^24,25,26 and for detecting genome-scale evolutionary patterns. Comparisons of the structure and sequence of these regions across different species within a genus can reveal important evolutionary phenomena such as gene transfer, deletion, or duplication. Recently, with the continuous application of high-throughput sequencing techniques, chloroplast DNA sequences have become readily available^13,27,28. However, there is no annotated chloroplast genome available for P. fraxinifolia, which hinders the understanding of the evolution of the chloroplast genome of this species from West Asia^14,15.

In this study, we aim to (1) assemble and annotate the chloroplast genome of the relict species P. fraxinifolia from Hyrcanian forests; (2) perform comparative genomics of the chloroplast genomes of six Pterocarya species; and (3) assess the systematic affinity of P. fraxinifolia using phylogenetic analysis of the assembled chloroplast genomes.

Leaf material for the P. fraxinifolia sample was collected from a wild population in Mazandaran, Iran (Fig. 1). The voucher samples were deposited at the Herbarium of the Nowshahr Botanical Garden (HNBG) under voucher number 12,876.

Genomic DNA was extracted using the CTAB method, and its quality and quantity were checked using a Qubit 2.0 and Agilent 2100 Bioanalyzer. Libraries were created and sequenced at Wuhan Benagen Tech Solutions Company Limited, Wuhan, China, using the DNBSEQ platform (paired-end 150 bp). SOAPnuke v1.3.0 was used to filter the raw data, yielding 20 GB of clean data²⁹.

Chloroplast genome assembly and annotation

Raw reads were filtered using Trimmomatic v0.39³⁰ with a quality cutoff of 15 in a 4-base sliding window; any reads that were less than 50 bp were removed, and the adapters were filtered out. The quality of the reads before and after trimming was tested using FASTQC v0.12.1. We used GetOrganelle³¹ v.1.7.7.0 for chloroplast genome assembly, with the embplant_pt database used as a reference and maximum extension rounds of 15 (-R). GetOrganelle produced two isomers of the whole chloroplast genome of P. fraxinifolia, and each genome had a distinct relative orientation for the small single-copy (SSC) region³². A Python script from GetOrganelle was used along with Bowtie2 v2.5.4³³ to determine the average read coverage throughout the chloroplast genome. GeSeq v2.03³⁴ was used for the initial chloroplast genome annotation of P. fraxinifolia, and the output from GeSeq was imported into Geneious Prime 2025.0.3 for an additional annotation check via the “Transfer Annotation” function. Chloroplot³⁵ was used to produce a circular representation of the plastome.

Comparative analyses of the Chloroplast genomes

Because the flanking inverted repeat (IR) regions of the chloroplast genome often vary among species, we used CPJSdraw³⁶ to compare the IR regions of the six species. We used CUSP from EMBOSS v6.6.0.0 to calculate relative synonymous codon usage (RSCU) for protein-coding genes of P. fraxinifolia. To identify simple sequence repetitions (SSR), we used a Perl script from the Microsatellite Identification tool (MISA)³⁷. The settings were adjusted to ten, five, and four repeats for mononucleotides, dinucleotides, and trinucleotides, respectively. Forward, reverse, palindrome, and complementary sequences with a minimum repeat length of eight bp and a maximum computed repeat of 50% were analyzed using REPuter³⁸. The complete chloroplast genome sequences of the six Pterocarya species were aligned with Fast Statistical Alignment v1.15.9³⁸ to perform the nucleotide diversity analysis. We used a Perl script (https://github.com/xul962464/perl-Pi-nucleotide-diversity) to estimate the nucleotide diversity (PI) with a sliding window analysis with a step size of 200 bp and a window length of 800 bp.

The selection pressure on chloroplast protein-coding genes (CDSs) was evaluated by aligning the nonredundant genes from six species using MAFFT v7.526³⁹. We ran ParaAT.pl v2.0⁴⁰ to compute synonymous substitution rates (Ks), nonsynonymous substitution rates (Ka), and Ka/Ks. Each CDS pair of one-to-one species combinations is used as a homolog with genetic code 11. We estimated Ka, Ks, and Ka/Ks among the six Pterocarya species with KaKs_Calculator v2.0⁴¹.

Phylogenetic analysis

We constructed a maximum likelihood (ML) phylogenetic tree to understand the relationships of Pterocarya species. Chloroplast genome sequences were acquired from GenBank for the other Pterocarya and related genera in the Juglandaceae family. The multiple sequence alignment contained a total of 21 taxa. We performed our phylogenetic analysis using the full chloroplast genome alignment, treating it as a standard coalescent gene⁴¹. The chloroplast genomes were aligned using Fast Statistical Alignment v1.15.9⁴² and then trimmed with trimAL v1.5⁴³ with the following settings: -automated1 -res overlap 0.7, -seqoverlap 65. To overcome the alignment issues, we also employed TAPER v1.0.047 with the -m N -a N parameters.

Using RAxML-NG v1.2.1⁴⁴, we constructed the GTR + G model and the ML tree with 500 bootstrap repetitions. The phylogenetic tree was rooted using Engelhardia roxburghiana Wall. as an outgroup. The tree was drawn using FigTree v1.4.4 (https://github.com/rambaut/figtree). To determine the genetic distance between the six Pterocarya species, the HKY85 model⁴⁵ was used, and a phylogenetic network was generated using the NeighborNet approach in SplitsTree CE v6.0.0⁴⁶.

Chloroplast genome assembly and annotation

The total numbers of raw and trimmed reads for P. fraxinifolia in this study were 143,190,876 and 141,927,817 base pairs (bp), respectively. The number of matched mapped pairs across the chloroplast genome was 393.42 ± 82.15 (Fig. S1). The complete chloroplast genome of P. fraxinifolia has a typical quadripartite structure that is 160,086 bp in length with a large single-copy region (LSC) of 89,582 bp, a small single-copy region (SSC) of 18,398 bp, and a pair of inverted repeat regions (IRs) of 26,053 bp (Fig. 2). A total of 148 genes were annotated in the chloroplast genome of P. fraxinifolia, including 103 protein-coding genes, 37 transfer RNA (tRNA) genes, and eight ribosomal RNA (rRNA) genes (Table 1 and Table S1). The GC content of the chloroplast genome was 36.17%. The annotated complete chloroplast genome of P. fraxinifolia was deposited in GenBank (accession number PV791734).

Comparative analyses of the Chloroplast genome and nucleotide diversity

According to a comparative analysis of the chloroplast genomes of Pterocarya species, the locations of eight genes in the chloroplast maps differed among species. The rps19 gene starts at position zero of the LSC region for P. fraxinifolia, but its position has shifted three times into the IRb region in the others. However, in other species of Pterocarya, a small portion of the genes were located in the IRb region. The ndhF gene in P. fraxinifolia, P. stenoptera, P. macroptera, and P. rhoifolia is located inside the SSC and is 2226 bp in length, whereas in P. tonkinensis and P. hupehensis, it spans 69 and 145 bp, respectively, into the IRb region (Fig. 3a).

The average nucleotide diversity (π) value was 0.001492, with a range of 0 to 0.00556 (Fig. 3B). The CDSs with the highest π values, which were greater than 0.0031, were ndhF, infA, ycf1, rps15, and matK. The ycf1 gene is found in the SSC area, whereas ndhF, infA, rps15, and matK are found in the LSC region. Nucleotide diversity decreased in both IR zones. Furthermore, 35 CDSs had a π value of zero among the six Pterocarya species, indicating that they were conserved (Table S1).

Repeated sequence analysis

The six Pterocarya chloroplast genomes have an average of 72.6 SSR loci (Fig. 4A), with P. rhoifolia having the most SSR loci (85) of the six species (Table S2). A thorough examination of the chloroplast genome of P. fraxinifolia revealed 68 microsatellites, comprising 63 mononucleotides, four dinucleotides, and one trinucleotide simple sequence repeat. The five types of sequence repeat motifs—forward, reverse, complementary, palindromic, and tandem—are summarized in Table S3 and Fig. 4B. The analysis also revealed that the number of repetitive sequences differed across the six Pterocarya chloroplast genomes. Approximately 96.82% of the mononucleotide repeats found in P. fraxinifolia were classified as A/T (61), and 3.18% (2 repeats) were classified as C/G. In contrast, approximately 88.2% of the repeats found in P. rhoifolia were classified as A/T (75), and 3.52% (3 repeats) were classified as C/G (Fig. 4C). Dinucleotide repeats (6) for P. rhoifolia and (4) for P. fraxinifolia were the next most prevalent type of SSR. This investigation revealed no repeats of tetranucleotides, pentanucleotides, or hexanucleotides.

Ka/Ks ratio and codon bias analysis

Strong purifying selection and functional limitations are indicated by the very low Ka/Ks ratios found in most CDS regions among Pterocarya species (Fig. 5A). With the exception of P. tonkinensis and P. stenoptera, the highest Ka/Ks ratio was detected in the chloroplast NADH dehydrogenase F (ndhF) gene. The GC contents for the first, second, and third codon locations were 45.30%, 38.25%, and 30.36%, respectively, whereas the overall coding GC content was 37.97%. The greatest frequencies were 42.361 for the ATT codon and 37.605 for the GAA codon. The only two codons with an RSCU value of 1 were tryptophan (TGG) and methionine (ATG) (Fig. 5B). Every codon ending in A or T had an RSCU value greater than 0.5.

Phylogenetic analysis

The aligned multiple sequence alignment for the phylogenetic analysis consisted of 158,422 bp across 21 accessions, with 0.21% gaps and 96.19% invariant sites. The phylogenetic tree revealed Pterocarya as a sister genus to Juglans L. with 100% bootstrap support (Fig. 6A). The ML phylogenetic tree confirmed the monophyly of the genus Pterocarya with 100% bootstrap support with two subclades. P. fraxinifolia is a sister to a monophyletic subclade that include P. tonkinensis and P. macroptera and a sister to another subclade that includes P. rhoifolia, P. stenoptera, and P. hupehensis. The network analysis of the six Pterocarya species revealed a topology similar to that of the ML tree, with P. tonkinensis clustering with P. macroptera and P. rhoifolia clustering with P. stenoptera and P. hupehensis, while P. fraxinifolia branched off independently. In this study, the efficiency of two barcode regions, matK and ycf1, in the phylogeny of the genus Pterocarya was evaluated (Fig. 6B and C). The results revealed that the phylogenetic tree based on the matK region was identical to the phylogenetic tree derived from the complete chloroplast genome sequence. Pairwise distance analysis using the HKY85 method revealed that P. fraxinifolia is distantly related to Asiatic Pterocarya species (Fig. S2). The genetic distances between P. macroptera and P. tonkinensis (0.000259) and between P. stenoptera and P. hupehensis (0.000526) were the lowest, whereas the genetic distances between P. fraxinifolia and P. hupehensis (0.002153) and between P. fraxinifolia and P. tonkinensis (0.001928) were more than eightfold greater (Table S4). In the MatK dataset, P. fraxinifolia had three unique character states that differentiated it from other species of Pterocarya (Table S5).

Chloroplast genomes are useful tools for studying the evolutionary relationships among species because of their preserved structure and uniparental inheritance (usually maternal in angiosperms^47,48. Considering mechanisms of plant evolution^49,50 and that the evolutionary history of chloroplasts is normally different from that of nuclear markers^51,52, the use of genetic information from chloroplasts could reflect how seed dispersal affects the genetic makeup of wild populations and species.

This study is the first to annotate the chloroplast genome of P. fraxinifolia and compare it to that of other species. We found that the positions of eight markers, namely, rps19, rpl2, ycf1 (IRa and IRb), ndhF, trnN, rpl2, and trnH, varied among the six Pterocarya chloroplast genomes. This implies that the expansion and contraction of the IR, LSC, and SSC areas are the primary sources of fluctuations in chloroplast genome size^53,54. Between 68 and 85 SSRs were found among the chloroplast genomes of the six Pterocarya species. While the number of poly(G)/(C) repeats was shown to be greater in other angiosperms, the number of poly(A)/(T) repeats was significantly greater in Pterocarya.

Five genes, ndhF, infA, ycf1, rps15 and matK, presented the greatest nucleotide variability (above 0.003). The matK and ycf1 genes have been suggested to function as barcode regions in plants⁵⁵. The matK gene encodes the maturase protein, which facilitates the splicing of group II introns in several chloroplast genes and is considered a core barcode for land plants^50,51. The ycf1 gene, which encodes the TIC214 protein that is essential for plant viability, is the second largest in the chloroplast genome and has recently been assessed for its DNA barcoding potential^50,51,52, showing higher variability than the existing chloroplast candidate barcodes (such as rbcL, matK and trnH-psbA). Therefore, the ycf1 gene might be potentially useful as a DNA barcode for the Pterocarya genus⁵⁶. With the exception of the matK region, none of the seven recommended barcode candidate genes in chloroplast genomes⁵⁰ have the potential for barcoding of the Pterocarya genus because of a lack of nucleotide variation. Surprisingly, the accuracy of the matK region in resolving the phylogeny of the genus Pterocarya was identical to that of the complete chloroplast genome. Therefore, the matK gene alone is sufficient for reconstructing the phylogenetic relationships within the genus Pterocarya, eliminating the need for the additional time and financial resources required for whole-chloroplast-genome sequencing.

The genus Pterocarya consists of six species and is closely related to Juglans in terms of pollen morphology, wood anatomy and molecular phylogenetics^8,9. Our phylogenetic results confirm the sister relationship of Pterocarya to Juglans. Two sections for Pterocarya have been proposed on the basis of the presence or absence of scales on the terminal buds^9,13,50. P. fraxinifolia, P. hupehensis, P. stenoptera, and P. tonkinensis belong to the section Pterocarya, while P. macroptera and P. rhoifolia belong to the section Platyptera. According to our chloroplast genome-based phylogeny, this suggested morphological classification is not supported, and the Caucasian wingnut (P. fraxinifolia) is in a distant subclade from the Chinese wingnut (P. stenoptera) and the Japanese wingnut (P. rhoifolia).

The pairwise genetic distance between the Caucasian wingnut and other Asiatic Pterocarya species is greater. This distance might reflect the prolonged isolation and considerable geographic distance between Caucasian wingnut and East Asian species. Recent divergence time analyses based on fossil calibrations estimated the age of P. fraxinifolia between 9.4 and 18.4 Ma from the Miocene period and suggested the westward dispersal of Pterocarya from East Asia⁸. Wingnut fruit structure could facilitate the dispersal of these species by wind and water⁵⁷. In this study, we collected P. fraxinifolia materials from its natural habitat in Hyrcanian forests. Our initial phylogenetic results revealed that the publicly available P. fraxinifolia in GenBank (NC046430) is not a P. fraxinifolia and is most likely a misidentified voucher that could be P. stenoptera (data not shown).

Toward conservation of P. fraxinifolia

P. fraxinifolia is classified as a vulnerable relict species on the IUCN Red List¹². Our phylogenetic tree, which was constructed on the basis of chloroplast genome analysis, indicates that this species is completely distinct from other species of the genus originating from China and Japan. This distinction might highlight the species’ unique evolutionary path and specialized ecological environments. Recent studies have shown that the potentially suitable ranges of P. fraxinifolia will increase under future climate scenarios^8,58, and the rapid loss of its habitat, combined with growing threats such as drought and the destruction of riparian ecosystems in Hyrcanian forests, will result in its conservation an urgent priority.