Four samples from China were sequenced in this study, and other genomes of the genus Torreya in East Asia were downloaded from GenBank database for comparison analysis (Table 1). Healthy young leaves were collected from each plant, wiped clean, numbered and placed in Ziplock bags. A large amount of color-changing silica gel was used for rapid drying.

The Hi-DNAsecure Plant Kit DP350 (Tiangen Biotech, Beijing, China) was used for DNA extraction following the manufacturer’s instructions (include manufacturer).

PE150 paired-end sequencing was conducted on an Illumina HiSeq X Ten genomic sequencer (Illumina Inc., San Diego, California, USA) at Majorbio Company (Shanghai, China). All sequencing depths were over 100$\times$. The subsequent analysis was based on clean reads with high quality.

MITObim v1.8 [30] was used to assemble the chloroplast genome with default parameters, and the final circular structure was formed manually. Annotations of these chloroplast genomes were performed by using online software GeSeq [31] and BLAST+ 2.9.0 [32] with an E-value of 1e-5. The genome map was illustrated by using OGDRAW [33]. Annotated chloroplast genome sequences were submitted to GenBank and raw data were uploaded to Sequence Read Archive (SRA) of NCBI and will be released upon publication. GenBank accession numbers of genomes used in this study are listed in Table 1.

Seven chloroplast genomes of the genus Torreya, including the four new genomes sequenced in this work, were compared by using mVISTA [34, 35] based on the Shuffle-LAGAN method [36].

The seven chloroplast genomes of Torreya species were also compared with BLAST Ring Image Generator (BRIG) v0.95 [37], and BLAST+ v2.9.0 was used with an E-value of 1e-5.

Single-copy orthologous genes were selected from the results analyzed by OrthoFinder v2.2.7 [38, 39]. We used MAFFT v7.45 for multiple sequence alignment with the L-INS-i strategy for more accuracy [40, 41]. Then, JModelTest v2.1.1.0 [42, 43] was used to find the best-fitting models for single-copy orthologous genes according to the Akaike information criterion (AIC), where the lowest value showed the best fit [44].

A phylogenetic tree based on ML with single-copy orthologous genes was built by using RAxML v8.2.4 [45] with the best-fitting model and 1000 bootstrap replicates. A BI-based phylogenetic tree was constructed with single-copy orthologous genes in MrBayes v3.2.7 [46, 47, 48], running for 2 $\times$ 10${}^6$ generations. The phylogenetic tree was illustrated by FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/) and TreeGraph v2.15.0-887 beta [49].

Results for the four sequenced genomes, including genome structure and GC content, are shown in Table 2. The genome sizes of T. parvifolia, T. nucifera, T. fargesii var. yunnanensis, and T. grandis var. jiulongshanensis are 136781 bp, 136955 bp, 136807 bp, and 137320 bp, with an average genome GC content of 35.49%, 35.46%, 35.49%, and 35.41%, respectively. Their large single-copy (LSC) regions are 97386 bp, 97516 bp, 97421 bp, and 114484 bp, respectively. Their small single-copy (SSC) regions are 38799 bp, 38841 bp, 38790 bp, and 22148 bp, respectively. The inverted repeat (IR) regions are 298 bp, 299 bp, 298 bp, and 344 bp, respectively. All of the genomes have similar small IR regions.

Chloroplast genome maps of T. parvifolia, T. nucifera, T. fargesii var. yunnanensis, and T. grandis var. jiulongshanensis are shown in Fig. 1. T. parvifolia, T. nucifera, and T. fargesii var. yunnanensis are very similar in gene number, order and names. However, T. grandis var. jiulongshanensis has a different gene number, order and name, especially gene order. In contrast to the other three taxa, T. grandis var. jiulongshanensis lacks the rps11 gene but has the clpP gene as a unique gene. The locations of some genes in the chloroplast genome maps differ from others. For example, in T. grandis var. jiulongshanensis, the atpA, atpF, atpH, atpI, rpoB, rpoC1, rpoC2, rps2, rps4, psaA, psaB, psbC and psbD genes are located in the LSC region, and the ndhF and ycf1 genes are located in the SSC region. However, in the other three species, the former group are located in the SSC region, and the latter group are located in the LSC region.

Fig. 1.

Chloroplast genome maps of four species of the genus Torreya in East Asia. (A) The chloroplast genome map of T. parvifolia. Genes outside the circle transcribe counterclockwise, while those inside transcribe clockwise. Genes with different functions are marked with different colors. (B) The chloroplast genome map of T. nucifera. Genes outside the circle transcribe counterclockwise, while those inside transcribe clockwise. Genes with different functions are marked with different colors. (C) The chloroplast genome map of T. fargesii var. yunnanensis. Genes outside the circle transcribe counterclockwise, while those inside transcribe clockwise. Genes with different functions are marked with different colors. (D) The chloroplast genome map of T. grandis var. jiulongshanensis. Genes outside the circle transcribe counterclockwise, while those inside transcribe clockwise. Genes with different functions are marked with different colors.

Seven chloroplast genomes of species in the genus Torreya, compared by mVISTA with the Shuffle-LAGAN method and annotated according to T. jackii MK249064, are shown in Fig. 2. Overall, the sequence similarity is high among these species, especially in exon regions. Regions with such divergence in exon regions include accD, ndhB, ndhF, psbA, psbJ, rpl2, rps3, rps16, rps18, ycf1, and ycf2. tRNA or rRNA genes with high levels of divergence include trnT-GGU, trnR-UCU, trnF-GAA, trnL-CAA, trnQ-UUG, trnN-GUU, trnM-CAU, rrn4.5 and rrn5. A comparison of the chloroplast genomes of seven species of the genus Torreya at the level of whole genome sequence was analyzed using BRIG, providing information for identifying the unique sequence of the chloroplast genome (Fig. 3). The sequence of T. jackii MK249064 was selected as a reference.

Fig. 2.

Comparison of chloroplast genomes of seven species of the genus Torreya by using mVISTA based on Shuffle-LAGAN method. The sequence of T. jackii MK249064 is selected as reference. The thick gray arrow at the top of the array indicates gene orientations. The dark-blue regions, light-blue regions, and pink regions represent exon, tRNA or rRNA genes, and conserved non-coding sequences (CNS), respectively.

Fig. 3.

Comparison of chloroplast genomes of seven species of the genus Torreya by using BRIG. The sequence of T. jackii MK249064 is selected as reference. Two of the innermost rings show GC content and GC skew from inside to outside, respectively.

Seventy-three single-copy orthologous genes were selected (Supplementary Table 1), and then a phylogenetic tree was built based on the ML and BI methods with these genes using the best-fitting model GTR + GAMMA + I. The ML and BI phylogenetic trees show the same relationships among the Torreya species (Fig. 4 ). T. parvifolia MN244711 is located closely to T. fargesii var. yunnanensis MN244712, and they are included in a clade with T. fargesii var. yunnanensis MK249061. T. fargesii NC_029398 and T. nucifera MK249060 form a single clade. T. parvifolia MN244711 and T. fargesii var. yunnanensis MN244712 are close to the clade containing T. fargesii NC_029398 and T. nucifera MK249060. T. grandis var. jiulongshanensis MN244714 is close to T. californica MK249062, and this clade is close to T. grandis NC_034806. The clade including T. grandis var. jiulongshanensis MN244714, T. californica MK249062 and T. grandis NC_034806 is close to our sample of T. nucifera MN244713. T. jackii KX902234 and T. jackii MK249064 are located in the most exterior part of the phylogenetic tree. The genus Torreya forms one clade in the phylogenetic tree. Based on phylogenetic trees, T. nucifera MK249060 is close to T. fargesii NC_029398, different from the results of Zhang et al. [1].

Fig. 4.

Phylogenetic tree shows relationship and supported values. Both ML and BI tree show the same relationship. Supported values of ML and BI tree are showed on branch and divide by sign of division, respectively.

Statistics of the genomes of four Torreya species showed that the IR regions were very short, with only 298–344 bp. Shrinkage or loss of the IR region is not rare in plants, and in a study of Pinus thunbergii, the IR region was found to have been reduced to 495 bp [50]. Loss of the IR region occurred in all cases of Pisum sativum L., Vicia faba L. [51], Glyptostrobus pensilis [52] and Cryptomeria japonica [53]. Pinaceae and non-Pinaceae conifers independently lost different copies of IRs, according to comparisons of the junctions near LSC regions and residual IR copies among gymnosperms [54], while the lack of an IR copy was considered a derived characteristic common to all conifers and “a single loss event defining the conifers as a monophyletic group” [55]. In addition, the four species shared similar GC contents (35.41%–35.49%). T. parvifolia, T. nucifera, and T. fargesii var. yunnanensis showed very similar chloroplast genome map structures, while T. grandis var. jiulongshanensis exhibited a very different order of genes. However, according to comparisons by both mVISTA and BRIG, there were no great differences among these species, especially in protein-coding regions. Therefore, this specificity could be due to genetic recombination.

In a 2006 publication, T. parvifolia was treated as a new species [13]. However, based on the phylogenetic analysis reported here, T. parvifolia MN244711 and T. fargesii var. yunnanensis MN244712 have a close relationship, with a high correlation and form one clade with T. fargesii var. yunnanensis MK249061. Therefore, we believe that the species T. parvifolia should be considered as T. fargesii var. yunnanensis or treated as a sister species according to the molecular evidence. T. nucifera MK249060 was collected from Lushan Botanical Garden, China, and the material of T. nucifera MN244713 used in this study was collected from Nanjing University, Nanking, Jiangsu Province (Jinling University introduced the species from Japan), which has been verified by historical evidence. In the recent study [1], T. nucifera MK249060 and T. fargesii NC_029398 formed a single cluster with a very high association based on the ML method, and this cluster was closest to T. fargesii var. yunnanensis MK249061, followed by T. grandis NC_034806. However, in this study, T. nucifera MN244713 was closest to T. grandis NC_034806 and T. grandis var. jiulongshanensis MN244714 and then to T. fargesii NC_029398 and T. fargesii var. yunnanensis MK249061. Therefore, the molecular biological evidence in this study supports the treatment of T. nucifera MK249060 as T. fargesii.

In a recent study by Kou et al. [15], phylogenetic analysis of the nuclear internal transcribed spacer (ITS) and combined sequences of chloroplast rbcL and rpl16 genes in the genus Torreya showed that T. grandis var. jiulongshanensis was close to T. jackii according to the combination sequences based on the ML method. However, it is insufficient evidence to suggest that T. grandis var. jiulongshanensis is a natural hybrid of T. jackii and T. grandis, basing on a few sequence fragments. In addition, the study comparing leaf variation among T. grandis var. jiulongshanensis, T. grandis, T. fargesii and T. jackii [16] proposed that T. grandis var. jiulongshanensis be treated as an independent species rather than a variety of T. grandis due to the morphological features of their leaves. The two studies above, examining gene fragments or just leaf variations, failed to offer hard evidence for related taxonomy. In our study, phylogenetic trees based on the whole genome-wide level single-copy orthologous genes showed consistent results by ML method and BI method, both illustrating that T. grandis var. jiulongshanensis MN244714 was close to T. grandis NC_034806. Therefore, according to our study, T. grandis var. jiulongshanensis should be treated as a variety of T. grandis rather than a natural hybrid between T. jackii and T. grandis or an independent species, supporting the original proposal [12].