We collected fresh young leaves of P. orientalis from an elite tree in a seed orchard in Jiaxian, Henan Province, China. DNA used for Illumina sequencing was isolated using a cetyl trimethyl ammonium bromide (CTAB)-based method (Doyle and Doyle, 1987). For the Oxford Nanopore Technologies (ONT) sequencing, we prepared DNA using a DNeasy Plant Mini Kit (QIAGEN). A short-insert 2 × 150 bp pair-end library was constructed following the manufacturer’s PCR-free protocol (Illumina) and sequenced on an Illumina NovaSeq 6000 platform. Library for ONT sequencing was prepared following the Nanopore 1D Genomic DNA by ligation (SQK-LSK109)–PromethION protocol, and sequenced on a PromethION platform. The short reads from Illumina were filtered using fastp v0.21.0 (Chen et al., 2018) with parameters of “-q 20 -l 36 --cut_right --n_base_limit 0 -c.” The long reads from Nanopore were corrected using NextDenovo v1.1.0 (https://github.com/Nextomics/NextDenovo). We generated ∼110 Gb Illumina clean short reads and ∼808 Gb (40.3 M reads with a N50 of 25 Kb) corrected ONT long reads.

We applied a long and short reads hybrid strategy to assemble the P. orientalis mitogenome. First, the clean short reads were de novo assembled using GetOrganelle v1.6.4 (Jin et al., 2020) with parameters of “-R 50 -k 67,87,107,127 -F embplant_mt --reduce-reads-for-coverage inf -w 127 -t 12 --max-reads 7.5E8”. The resulting unitig graph was manually edited using Bandage v0.8.1 (Wick et al., 2015) to remove chloroplast- and nuclear-derived unitig nodes. The corrected long reads were then mapped to the graph using the minimap2 v2.24-r1122 (Li, 2018). Finally, repeats on the graph were resolved by aligning with the mapped long reads.

The protein-coding genes, non-coding rRNA, and tRNA genes were annotated using the Organelle Genome Annotation Pipeline (OGAP, https://github.com/zhangrengang/OGAP) with parameters of “-taxon Acrogymnospermae Bryophyta Liliopsida commelinids fabids malvids -mt -trn_struct”. More specifically, protein-coding genes were identified with Exonerate v2.2.0 (Slater and Birney, 2005) and AUGUSTUS v3.3.3 (Stanke et al., 2008), tRNA genes with tRNAscan-SE v2.0.5 (Lowe and Eddy, 1997), and rRNA genes with BLAT v36x2 (Kent, 2002) within the OGAP.

RNA editing sites in protein-coding genes of each species were predicted using the PREP-Mt (Mower, 2005) with a cutoff of 0.2. To validate the accuracy of the predictions, we searched for RNAseq data of P. orientalis (Supplementary Table S1) and mapped them to the protein-coding genes with bwa mem v0.7.17-r1188 (Li, 2013) to determine the edited sites using REDItoolDenovo v1.3 (Lo Giudice et al., 2020) with minimum coverage of 10 and mapping quality score of 30.

To estimate the repetitive sequences in mitogenomes, we annotated tandem repeats, TEs, and dispersed repeats in P. orientalis and 19 other taxa representing the six major lineages of seed plants included in this study (Table 1). Tandem repeats were identified using Tandem Repeats Finder v4.10.0 with the parameters of “2 7 7 80 10 50 800 -d” (Benson, 1999). For TEs, we first used the de novo TE family identification and modeling package RepeatModeler v1.0.10 (Smit and Hubley, 2008) to train a TE database from each mitogenome. The RepeatModeler derived library was then combined with the database “RepeatMasker.lib” of RepeatMasker v4.07 (Smit et al., 2013) and the repeat library of Picea abies v1.0 assembly (Nystedt et al., 2013) to generate a joint repeats library. Finally, TEs in each genome were identified by RepeatMasker using this reference database.

Dispersed repeats were identified with a BLASTN search using a word size of seven and an e-value of 1E-6 as in a previous study (Guo et al., 2016). Overlapping regions of dispersed repeats were counted only once when calculating their total size. The dispersed repeats were grouped into small-, medium-, and large-size classes with sequence lengths of <100 bp, 100–1,000 bp, and ≥1,000 bp, respectively. Divergence among repeats was estimated using Kimura two-parameter distance (K2P) (Kimura, 1980) for pairwise repeat sequences.

Based on gene annotation results, we found 8–12 genes were absent in mitogenomes of Cupressophytes and Gnetophytes (see Section 3.3). To investigate the possibility of their transfer to nuclei, we first tried to recover these genes from transcriptomes. We searched for RNAseq data of Cupressophytes and Gnetophytes in NCBI (https://www.ncbi.nlm.nih.gov/), and de novo assembled transcriptomes for Cupressus sempervirens, Hesperocyparis glabra, Taxus cuspidata, Podocarpus macrophyllus, Gnetum gnemon, and Welwitschia mirabilis (Supplementary Table S2). The Illumina reads were filtered using fastp v0.21.0 (Chen et al., 2018) with the parameters of “-q 20 -l 36 --cut_right --n_base_limit 0 -c,” generating 6.3–14.7 Gb clean reads per species. We then de novo assembled the transcriptomes using Trinity v2.8.5 with default parameters (Grabherr et al., 2011). The transcriptome of P. orientalis is from our early study (Hu et al., 2016). For Araucaria heterophylla, we only found 212 ESTs/cDNAs in TreeGenes database (Wegrzyn et al., 2008) (Supplementary Table S2). We searched for intact mitochondrial genes in these transcriptomes and ESTs/cDNAs using OGAP with the parameters of “-taxon Acrogymnospermae Bryophyta Liliopsida commelinids fabids malvids -mt -trn_struct -trans”. Potential chloroplast genes were filtered out after BLAST against NCBI’s NT database (E-value < 1E-40, sequence coverage >99%, and sequence identity >99%).

Second, to confirm that the transcriptome-recovered mitochondrial genes (mitogenes) were in nuclear genomes, we used TBLASTN to search for matches of these genes in reference genomes of Cupressophytes and Gnetophytes, including Sequoiadendron giganteum (Scott et al., 2020), Sequoia sempervirens (Neale et al., 2021), Taxus chinensis (Xiong et al., 2021), Taxus wallichiana (Cheng et al., 2021), Gnetum montanum (Wan et al., 2018), and Welwitschia mirabilis (Wan et al., 2021) (Supplementary Table S3).

Third, we used mapping depth of reads of resequencing data to verify the genome origins of the transcriptome-recovered mitogenes. In a plant cell, the abundance of nuclear and organelle genomes can differ by thousands fold with the copy numbers of plastid the highest, mitochondrial intermediate and nuclear the lowest (Petit and Vendramin, 2007). Based on this expectation, we investigated the depth of read coverage of each recovered genes. The mitochondrial genes retrieved from mitogenome and transcriptomes were used as reference. We also include the low-copy nuclear gene LEAFY and chloroplast gene matK as spike-in reference. We downloaded whole-genome resequencing data of Cupressus sempervirens, Hesperocyparis glabra, Taxus cuspidata, Podocarpus macrophyllus, Gnetum gnemon, and Welwitschia mirabilis from Sequence Read Archive (SRA) of NCBI (https://www.ncbi.nlm.nih.gov/sra) (Supplementary Table S4), the resequencing data of P. orientalis is from this study. The Illumina reads were filtered using fastp v0.21.0 (Chen et al., 2018) with the same parameters as in Section 2.1, generating 7.2–31.4 Gb clean reads per species. The clean reads were mapped to the reference genes using bwa mem v0.7.17-r1188 with default parameters. The resulting sequence alignment maps (SAMs) were sorted with SAMtools v1.14 (https://github.com/samtools/samtools), reads with alignment score (AS) <25 were removed by using ngsutilsj (https://github.com/compgen-io/ngsutilsj). We used SAMtools to compute the depth of the aligned reads at each position of the gene.

To identify synteny blocks among the mitogenomes of the major gymnosperm lineages, we performed pairwise genome alignment using Mauve v2.4.0 with default parameters (Darling et al., 2010). To calculate the degree of shared DNA between species, each pair of mitogenomes were aligned using BLASTN with a word size of seven and an E-value of 1E-6.

We used 28 protein-coding genes shared among all 20 mitogenomes included in this study to reconstruct a phylogenetic tree. Edited protein sequences were individually aligned with MAFFT v7.453 with default parameters and converted into the corresponding edited codon sequence alignments. Poorly aligned regions were trimmed using BMGE v1.1 (Criscuolo and Gribaldo, 2010) with parameters of “-t CODON -g 0.5.” The trimmed alignments of each gene were concatenated into a single alignment of 26,844 bp. This alignment was used to build a maximum-likelihood (ML) phylogenetic tree using IQ-TREE v1.6.12 (Nguyen et al., 2015), using the best-fit model TVM + F + R3 selected by ModelFinder (Kalyaanamoorthy et al., 2017) and with 1,000 replications of ultrafast bootstrap and Shimodaira-Hasegawa-like approximate likelihood-ratio (SH-aLRT) test. The branch lengths in units of synonymous (dS) and nonsynonymous (dN) substitution rates were estimated under the free-ratio branch model using codeml in PAML v4.9j package (Yang, 2007).

To estimate the divergence time between lineages, we run the MCMCTree in PAML v4.9j (Yang, 2007) for 2,000,000 iterations and 100,000 iterations for a burn-in. To calibrate the tree, we used the following time points of diversification from the TimeTree web (http://www.timetree.org/): 289–337 MYA between gymnosperms and angiosperms, 282–324 MYA between Cycads and other gymnosperm clades, 159–246 MYA between Cycas taitungensis and the other Cycads species, 134–197 MYA between Abies sibirica and other species within Pinaceae, 7.6 to 30.7 MYA between Picea abies and Picea glauca, 59–87 MYA between Pinus taeda and Pinus lambertiana, 99–123 MYA between Gnetum gnemon and Welwitschia mirabilis, and 213–272 MYA between the clades consisting of Podocarpus macrophyllus and Araucaria heterophylla and the other species of Cupressophytes.

We assembled the mitogenome of P. orientalis into two circular-mapping chromosomes representing a master ring and a sub-genomic ring with a total size of ∼2.6 Mb (Figure 1A, Table 2). This genome lacks large structural variants but is rich in dispersed repeats (Figures 1B,C). The genome size is similar to C. sempervirens (2.7 Mb) and H. glabra (2.0 Mb) in the family of Cupressaceae, but more than 5-fold larger than T. cuspidata and ∼2–4 fold smaller than Picea and Larix in Pinaceae (Table 1). Across the Cupressophytes and Pinaceae, the mitogenomes vary 6–10 fold in size, in contrast to the relatively stable and small genomes in Cycads and Gnetophytes.

The mitogenome of P. orientalis contains 32 distinct protein-coding genes (Figure 1A, Table 2), seven tRNA genes (Asp-trnD-GUC, Glu-trnE-UUC, Ile-trnI-CAU, Met-trnM-CAU, Met-trnfM-CAU, Trp-trnW-CCA, and Tyr-trnY-GUA), and three rRNA genes. Among them, ccmFC and cox1 have four copies, rrn5, rrn18, and rrn26 have three copies, nad1, trnD-GUC, trnE-UUC, trnM-CAU, and trnfM-CAU have two copies. We found five cis-spliced introns and 10 trans-spliced introns in the protein-coding genes (Supplementary Table S5). Both the number of protein-coding genes (Figure 2B) and the number of cis- and trans-spliced introns identified in P. orientalis are similar to that of Cupressophytes taxa reported previously (Guo et al., 2020). We also detected seven complete chloroplast genes (although they may not be functional), atpA, ndhA, ndhB, ndhC, ndhD, ndhE, and ndhH, and two chloroplast pseudogenes, petB and petD, in the mitogenome of P. orientalis (Table 2).

RNA editing remodels variability in higher plant mitogenomes, and could bias phylogenetic inferences if not accounted for (Szmidt et al., 2001; Mower, 2005). We predicted the RNA editing sites in protein coding genes of all 20 species and found that Ginkgo (1,300 sites) and Gnetophytes (225–336 sites) had the highest and lowest C-to-U editing sites across all protein-coding regions, respectively, and the other species spanned a range of 520–1,220 sites (Figure 2C). This variation in editing frequency reflects a lineage-specific pattern of editing, for example 25% (254) and 44–45% (102–147) of the editing sites in Ginkgo and Gnetophytes were unique (Figure 2C). In contrast, C. sempervirens and H. glabra of Cupressophytes each had only 6% (43–47) unique sites (Figure 2C).

To assess the accuracy of the predicted editing sites, we identified the edited sites from RNAseq data in P. orientalis. A total of 1,068 C-to-U edit sites were detected in 32 genes, of which 745 sites were among the 1,007 predicted sites (Table 2, Supplementary Table S6). This comparison suggests that the prediction method provides a reasonably reliable estimate of RNA editing (74% in P. orientalis), consistent with the results from a previous study (Guo et al., 2016).

We inferred the phylogeny of the major gymnosperm lineages using 28 conserved protein-coding genes common to all 20 species included in this study. The concatenated sequence alignment consisted 26,844 sites. On this phylogenetic tree, Cupressophytes and Gnetophytes are sister clades with a strong bootstrap support of 92% (Figure 2A). Within Cupressophytes, P. orientalis is sister to the group consisting of C. sempervirens and H. glabra (Figure 2A), and diverged from this group 87 MYA in the Paleogene. The divergence of the Cupressaceae family is dated to 145 MYA (Figure 2A), which is similar to the 159 MYA estimate from nuclear data (De La Torre et al., 2017).

To evaluate the substitution rates in mitochondrial genes, we estimated synonymous (dS) and nonsynonymous (dN) substitution rates along the branches of the inferred phylogeny. We found that dS and dN varied up to 46- and 234-fold among species, respectively. This wide range is mainly caused by Cupressophytes and Gnetophytes with their substantially longer branches and larger variation in dS and dN among species (Figure 2D; Table 3) compared to the other clades, all of which have low substitution rates. The overall dS of Cupressophytes and Gnetophytes is more than 5-fold higher than Pinaceae. The same overall pattern is seen with dN, although it is about 20–70% of the dS values, reflecting selective constraints on nonsynonymous mutations (Figure 2D; Table 3). Within Cupressophytes, dS varied by 46-fold and dN varied by 25-fold among lineages, with P. orientalis at the low end for both values (dS = 0.0030, dN = 0.0015; Table 3).

Moreover, we calculated the absolute substitution rates RS and RN per branch by adjusting the dS and dN by dividing the divergence time of the branch. Again, RS and RN varied up to 118- and 171-fold among species, respectively, with the highest RS and RN in Cupressophytes and Gnetophytes, and the lowest in Pinaceae. Within Cupressophytes, RS and RN varied by 40-fold and 70-fold, with P. orientalis having the lowest rates (Table 3). This wide range of substitution rates indicates that sequence evolution is heterogeneous within and among plant lineages, especially in Cupressophytes.

Early studies suggest that the mitogenome of the common ancestor of seed plants probably contained 41 protein-coding genes (Richardson et al., 2013; Guo et al., 2016). We found all 41 genes in Cycads, Ginkgo and Pinaceae, but 8–12 genes were absent in Cupressophytes and Gnetophytes (white and orange boxes in Figure 2B). All species in Cupressophytes have lost the same set of eight ribosomal protein genes, while Cupressaceae species and T. cuspidata lost an additional sdh3 gene (Figure 2B). In Gnetophytes, both G. gnemon and W. mirabilis have lost 11 genes, including 10 ribosomal protein genes and the sdh3 gene; their rpl16 appeared as degraded pseudogene (Figure 2B).

The reduction in mitochondrial protein-coding genes in Cupressophytes and Gnetophytes has been suggested to be a result of gene transfer to the nucleus prior to their loss from the mitogenome (Guo et al., 2016; Guo et al., 2020). To verify whether the lost genes were transferred to the nucleus, we first searched for them in the transcriptomes of Cupressophytes and Gnetophytes (Supplementary Table S2). We recovered 4–7 intact lost mitochondrial genes in transcriptomes of each respective species, but none in A. heterophylla (Figure 2B). The recovery of few intact mitochondrial genes in A. heterophylla is likely to be a consequence of the limited numbers of ESTs or cDNA sequences in the database, which would hamper detection of these genes in its nuclear genome. We then searched for these recovered genes in their own or related species’ reference genomes and found all of them with high identity (89–100%) and coverage (85–100%) scores (Supplementary Table S3). There were a few exceptions, e.g., rps11 and sdh3 showed lower identity (52–75%) but high coverage (92–100%) scores, likely due to different substitution rates between genes and species. The BLAST hits of Podocarpus macrophyllus genes gave moderate identity scores (48–81%) likely because of the lack of closely related reference. We further analyzed the mapping depth of genome resequencing reads to the mitogenes and transcriptome-recovered genes in each Cupressophytes and Gnetophytes species. In P. orientalis, the median depth of read coverage of the 32 genes in mitogenome ranged 275–1,169×, while the nuclear gene LEAFY had 15× and the chloroplast matK 14,395× (Supplementary Figure S1). The six recovered mitogenes had median depth of read coverage 11–17×, similar to the LEAFY, supporting their location in nucleus. Using the same approach, we confirmed that 4–7 of the lost genes in the seven other species of Cupressophytes and Gnetophytes were transferred to the nucleus (orange boxes in Figure 2B). Our results illustrate the differential abundance of nuclear, mitochondrial and plastid genomes in plant cells, and demonstrate the power of depth of read coverage as a quantitative tool for verification of genome origin of a gene.

Rapid structural evolution is thought to be characteristic of plant mitogenomes (Palmer and Herbon, 1988; Cole et al., 2018). To understand the rate of sequence turnover and structural rearrangement in gymnosperms, we examined the degree of synteny and shared DNA among the mitogenomes of the main gymnosperms lineages. We found that syntenic blocks and shared DNA between genomes reduced rapidly with phylogenetic distance (Figures 3A–D; Table 4). For example, P. orientalis and C. sempervirens diversified for 87 million years (MY) shared ∼60% DNA sequence, but this value dropped sharply to only ∼10% between P. macrophyllus and A. heterophylla with a divergence time of 126 MY (Figures 3A,C; Table 4). We examined the patterns of shared DNA between 16 pairs of gymnosperm species, and found a strong negative correlation (r = 0.53, p = 0.03) with time of divergence (Figure 3F). However, as an outlier, exceptionally low DNA turnover was found between G. biloba and C. taitungensis (Table 4), similar to the report by Guo et al. (2016). In addition, high proportions of shared DNA (58–64%) were retained between C. taitungensis and Ceratozamia hildae/Zamia integrifolia with the divergence time of 194 MY (Figure 3E; Table 4). These results illustrate a variable rate of mitogenome turnover in gymnosperms.

Repetitive sequences, particularly TEs, are important components of plant nuclear genomes and contribute to genome size variation and functional regulation (Mehrotra and Goyal, 2014). In contrast, TEs and other repetitive elements in plant mitogenomes are poorly characterized. In this study, we annotated TEs, tandem repeats, and dispersed repeats in the 20 mitogenomes, and found that they have a strong correlation with the genome size and are highly variable among species (Figure 4; Table 5). The highest proportion of TE is in A. sibirica (19.26%), followed by Z. integrifolia (15.57%), and the lowest in W. mirabilis (1.06%). TE contents vary widely even in the same clade, e.g. from 1.58% (P. macrophyllus) to 6.80% (H. glabra) in Cupressophytes, from 5.29% (P. lambertiana) to 19.26% (A. sibirica) in Pinaceae, and from 2.52% (C. taitungensis) to 15.57% (Z. integrifolia) in Cycads. Among TEs, LTRs (long terminal repeats) are the most common group followed by LINE (long interspersed nuclear elements) in most species (Figure 4A). For the three Cupressaceae species and the two Pinus species, the most abundant element is the “Unknown” category (Figure 4A). The “Unknown” category of TEs could be due to accumulation of various mutations over long evolutionary time that have degraded the TEs into sequences that cannot be properly classified.

In addition to TEs, dispersed repeats, including inverted and direct orientation duplicated sequences, are prevalent and made up to 20–30% of the mitogenomes, with a few exceptions (Table 5). The highest proportion is found in Cycads (22–32%) and A. sibirica in Pinaceae (32%). The two species in Gnetophytes are highly variable in this regard with 6% in W. mirabilis but 27% in G. gnemon; this value was stable (21–22%) among the three taxa in Cupressaceae. Tandem repeats are a smaller component compared to dispersed repeats and TEs, and were observed at the level of 1–5% in most of the mitogenomes (Table 5). The overall repetitive sequences composed 20–38% of the analyzed mitogenomes with a few exceptions of 9–17%, and show a strong correlation with genome size (r = 0.99, p < 0.0001; Figure 4E; Table 5).

To understand the dynamics of dispersed repeats in mitogenomes, we divided the repeats into small-, medium-, and large-size classes with lengths of <100 bp, 100–1,000 bp, and ≥1,000 bp, respectively. We found that small-size repeats are the most abundant class in the 20 mitogenomes, followed by the medium-size repeats, while large-size repeats are rare (Table 5). The median lengths for the three classes are 52 bp, 136 bp, and 1,517 bp, respectively. The number of small-size repeats are 2–20 fold more common than medium-size repeats; G. biloba was extreme having 112-fold more small-size than medium-size repeats (Table 5). In A. sibirica, the two classes are balanced. Large-size dispersed repeats are limited, with 50% of the species containing fewer than 15 large repeats (Table 5). On the other hand, three species A. sibirica, H. glabra and L. sibirica, had more than 100 large repeats (Table 5). Thus, the density plots for small-size and medium-size repeats are similar among species, while it varied widely for the large-size repeats (Figure 5A). The GC contents are relatively conserved for the three classes of dispersed repeats (Figure 5B).

To understand the tempo of evolution in dispersed repeats, we analyzed the distribution of K2P distances among the repeats in the 20 mitogenomes. Small values of K2P represent recent duplications while large values represent more ancient duplication events. We divided the K2P values into three distance classes: ∼0 as recent duplications, >0 and <0.2 as intermediate events, and >0.2 as ancient insertions. A K2P distribution with a dominant narrow peak indicates rapid accumulation of repeats in a short period, otherwise a flattened peak indicates a slow generation of repeats over a prolonged period. As shown in Figure 5C, the K2P distributions of the small-size class is similar among all species with a dominant recent peak and lack of ancient insertions, suggesting they have been accumulated recently over a short time span. In contrast, the K2P distribution of the medium-size and large-size repeats are visibly different among species, often with two or three peaks spanning recent to ancient classes, suggesting multiple and distant expansion events.

Ancient dispersed repeats are easy to identify based on K2P distances, but they are not prevalent in the 20 mitogenomes. This may suggest that the persistance of ancient large-size repeats is governed by the rates of accumulation and elimination. Large repeats can be eliminated or broken down by repeat-mediated recombination, resulting in an increase in small- or medium-size repeats that appear to be more recent expansions (Maréchal and Brisson, 2010; Sullivan et al., 2020).

Among the available gymnosperm mitogenomes, genome sizes vary markedly by 30-fold (Table 1). This large variation is hypothesized to be due to lineage- or species-specific accumulation of repetitive sequences, horizontal transfers and intergenomic transfers (Goremykin et al., 2012; Rice et al., 2013; Guo et al., 2017). It has been reported that 29–100% of sequences in gymnosperm mitogenomes are similar to the nuclear genome sequences (Sullivan et al., 2020). In P. orientalis we found 76% undefined sequences (Table 2), which could in part have been derived from the nucleus by intracellular transfer.

Another source of mitogenome size variation is repetitive sequences. Estimations of repetitive contents can be affected by filter and search criteria, e.g. setting the length to 50 bp or 100 bp could underestimate the full size range of repeats. In this study, we searched for dispersed repeats in all size classes. We found that dispersed repeats occupied 6–32% and TEs 1–19% in the sampled mitogenomes. The wide range of dispersed repeats content quantified in this study mirror the results of 1–49% found in 82 angiosperms (Dong et al., 2018). All categories of repetitive sequences are strongly correlated with genome size, giving an overall relationship of r = 0.99. Our results illustrate that repetitive sequences, especially TEs and dispersed repeats, are major components of mitogenome expansion.

Although a major component of mitogenomes, the size distributions of repetitive sequences and the dynamics of their accumulation have rarely been investigated in gymnosperms. We analyzed the frequency and sequence diversity of different size classes of repeats, and found small-size repeats are most common while large-size repeats are rare and vary widely among species. K2P distances among repeats indicate that most of the repeats have accumulated by recent expansions. Large-size repeats usually tend to recombine more frequently than small repeats (Skippington et al., 2015; Guo et al., 2016), which may eliminate large repeats over time. However, in P. abies, most recombinogenic repeats are small- to medium-size, which deviates from the general expectation that, in the absence of other factors, recombination should scale positively with repeat length (Sullivan et al., 2020). We are of the opinion that genome- or region-specific recombination rates could lead to unique sequence turnover rates resulting in the heterogeneous distribution of large repeats in mitogenomes.

The plant mitogenome is characterized by rapid structural rearrangement and sequence turnover, driven by a high recombination activity, making mtDNA highly variable even between closely related species (Guo et al., 2016; Cole et al., 2018; Sullivan et al., 2020). In line with this view, extensive rearrangement is found between two spruce species, P. abies and P. glauca (Sullivan et al., 2020). On the other hand, the mitogenomes of G. biloba and C. taitungensis have retained substantial amounts of shared DNA despite their divergence time >300 MY (Table 4; Guo et al., 2016). Even with these extreme cases, we found that the synteny blocks and shared DNA between mitogenomes generally decayed with divergence time, with a negative correlation r = 0.53 (p = 0.03) among the gymnosperms (Figure 3F). Although not exhaustive, our analysis provides a quantitative estimation of the rate of sequence turnover in gymnosperm mitogenomes. Previous studies show that rearrangement in mitogenomes occurs largely through repeat-mediated recombination (Cole et al., 2018; Xia et al., 2020). However, the relationship between inter- and intra-molecule recombination and rearrangement rates among species is unclear (Cole et al., 2018). Recombination dynamics in plant mitogenomes are poorly understood to date and deserve further scrutiny as it drives sequence and structural evolution of plant mitogenomes, and affects genome integrity and biological function.

The phylogenetic placement of Gnetophytes differs among studies. Five conflicting hypotheses were supported by different marker systems and include: the Gnetophytes-other seed plant hypothesis, the Gnetophytes-other gymnosperm hypothesis, the Gnetifer (Gnetophytes + Conifers) hypothesis, the Gnepine (Gnetophytes + Pinaceae) hypothesis, and the GneCup (Gnetophytes + Cupressophytes) hypothesis (Ran et al., 2018). Yet another suggestion is that Gnetophytes is sister to or within the conifers (Wan et al., 2018). The Gnepine hypothesis is favored by most studies (Ran et al., 2018), including a recent phylogeny based on mitochondrial protein-coding genes (Kan et al., 2021). On the other hand, plastid DNA strongly supports the GneCup hypothesis (Ruhfel et al., 2014; Gitzendanner et al., 2018; Leebens-Mack et al., 2019). This incongruence between datasets might be due to differences in the substitution rates of mitochondrial, plastid and nuclear genes, as well as the amount of informative sites included in each study.

In this study, we inferred the relationships among major gymnosperm lineages using 28 conserved protein-coding genes. The recovered phylogeny supports the GneCup hypothesis with a strong bootstrap support. We are more in favor of GneCup hypothesis because of the shared gene transfer pattern in these two clades. Over evolutionary time, some segments of mtDNA have diverged so much that groups of genes have been lost. The common ancestor of seed plants is inferred to contain 41 protein-coding genes (Guo et al., 2016; Mower, 2020), and Cycads, Ginkgo, and Pinaceae maintained this same set. The mitogenomes of Cupressophytes and Gnetophytes stand out from other gymnosperms in that both have lost 8–11 ribosomal protein-coding genes and the sdh3 gene, with the exception of Ephedra przewalskii (Gnetophytes) which has lost 19 genes (Guo et al., 2020; Kan et al., 2021). To confirm whether these lost genes have been transferred to nucleus, we first annotated them in transcriptomes, then identified their homologs in nuclear genomes, and further examined their depth of mapped reads using genome resequencing data. All three lines of evidence support the transfer of 4–7 mitogenes to nucleus. Our findings are in good agreement with the recent studies by Kan et al. (2020) and Kan et al. (2021), in which they report finding 6–8 intact or partial homologs of lost genes in the nuclear genomes. In sharp contrast, the three other lineages of gymnosperms, including Pinaceae, Ginkgo, and Cycads, share the same set of genes as the common ancestor more than 300 MYA (Guo et al., 2020). The unrecovered missing mitogenes could either be truly lost or not present in the transcriptomes we analyzed due to tissue- and/or time-specific expression, or incomplete annotation from transcriptomes, which identified the genes as partial homologs instead of intact genes.

Other striking characteristics of Cupressophyte and Gnetophyte mtDNA are their accelerated and hugely variable substitution rates across the 28 conserved protein-coding genes relative to other gymnosperms. The synonymous substitution rate dS in Cupressophytes and Gnetophytes averaged 0.0641 and 0.0536, respectively. After scaling up by divergence time, the absolute synonymous substitution rates RS were similar in Cupressophytes and Gnetophytes with 2.5 × 10⁻¹⁰ and 2.1 × 10⁻¹⁰ site/year, respectively, and are ∼2–7 fold higher than Pinaceae, Ginkgo, and Cycads. Previous work similarly inferred exceptionally low substitution rates in Pinaceae (Pinus and Picea) and higher substitution rates in Gnetophytes species (Guo et al., 2016; Sullivan et al., 2020; Kan et al., 2021). However, the extraordinary variation in RS in Cupressophytes and Gnetophytes is revealed only by denser sampling of taxa from each clade. Within Cupressophytes, RS varied 46-fold among taxa. Our results add to the recognition that while most plant mitogenomes show very low synonymous substitution rates, this rate of divergence is also exceptionally variable for reasons still unclear (Mower et al., 2007). Another unexpected finding is that within Gnetophytes, G. gnemon and W. mirabilis showed distinct species-specific RNA editing patters as 87–91% of their editing sites differ from each other. Such a striking diversification in RNA editing is not observed within other clades. Future studies with broader taxa sampling are required to understand the mechanisms driving these patterns, including gene loss and transfer, TEs, recombination and selection.