Revisiting the Plastid Phylogenomics of Pinaceae with Two Complete Plastomes of Pseudolarix and Tsuga.

Phylogeny of the ten Pinaceous genera has long been contentious. Plastid genomes (plastomes) provide an opportunity to resolve this problem because they contain rich evolutionary information. To comprehend the plastid phylogenomics of all ten Pinaceous genera, we sequenced the plastomes of two previously unavailable genera, Pseudolarix amabilis (122,234 bp) and Tsuga chinensis (120,859 bp). Both plastomes share similar gene repertoire and order. Here for the first time we report a unique insertion of tandem repeats in accD of T. chinensis From the 65 plastid protein-coding genes common to all Pinaceous genera, we re-examined the phylogenetic relationship among all Pinaceous genera. Our two phylogenetic trees are congruent in an identical tree topology, with the five genera of the Abietoideae subfamily constituting a monophyletic clade separate from the other three subfamilies: Pinoideae, Piceoideae, and Laricoideae. The five genera of Abietoideae were grouped into two sister clades consisting of (1) Cedrus alone and (2) two sister subclades of Pseudolarix-Tsuga and Abies-Keteleeria, with the former uniquely losing the gene psaM and the latter specifically excluding the 3 psbA from the residual inverted repeat.

Introduction

Pinaceae, the largest family of conifers, comprises more than 230 species in 10 genera—Abies Mill., Cathaya Chun & Kuang, Cedrus Trew, Keteleeria Carrière, Larix Mill., Picea A. Dietr., Pinus L., Pseudotsuga Carrière, Pseudolarix Gordon, and Tsuga (Endl.) Carrière. The family is an important resource for timber, pulp, essential oils, and other forest products. The Pinaceae are exclusively distributed in the northern hemisphere, except for one species, Pinus merkusii Jungh. & de Vriese, whose habitat crosses the equator in Sumatra (Thieret 1993).

The plastid genomes (plastomes) of photosynthetic seed plants are typically small (∼150 kb) with a quadripartite structure containing two inverted repeats (IRA and IRB, ∼20 to 30 kb each), which separate the large and small single copy regions (LSC and SSC) (Jansen and Ruhlman 2012). However, the plastomes of Pinaceous species only range from 107 to 120 kb (Lin et al. 2010) because of their highly reduced IRs (Wu et al. 2007; Wu, Wang, et al. 2011). In addition, Wu, Lin, et al. (2011) reported four distinct plastomic organizations among Pinaceous genera. The diversity of Pinaceous plastomic forms was proposed to be associated with intraplastomic homologous recombination, which is mainly triggered by two types of Pinaceae-specific IR, type 1 and 3 repeats (Wu, Lin, et al. 2011).

The plastome has served as a practical resource to resolve many questions in evolutionary studies, especially in green plant phylogeny (Ruhfel et al. 2014). Previously, Pinaceae plastid phylogenomic study (Lin et al. 2010) evaluated the phylogenetic relationships among eight of the ten Pinaceous genera. However, the remaining two, Pseudolarix and Tsuga, were not included in the study because of the unavailability of samples. Pseudolarix is a monotypic genus restricted to hills and plains along the Yangtze River valley in southeast China (LePage and Basinger 1995), whereas Tsuga contains nine recognized species in East Asia and North America (Havill et al. 2008).

In this study, we determined the complete plastomes of Pseudolarix amabilis and T. chinensis. Comparative plastomic analyses across the ten Pinaceous genera revealed that the accD of Tsuga is expanded with tandem repeats of PD/H amino acids. Our plastid phylogenomic results indicate that the five genera of Abietoideae constitute a monophyletic clade with Cedrus as sister to the clade of the other four genera, including Abies, Keteleeria, Pseudolarix, and Tsuga; and Pseudolarix and Tsuga form a subclade as a sister to the Abies-Keteleeria subclade.

Results and Discussion

Plastomic Features of P. amabilis and Tsuga chinensis

The plastomes of P. amabilis (LC095867) and T. chinensis (LC095866) are circular molecules of 122,234 and 120,859 bp, respectively (supplementary fig. S1, Supplementary Material online). Like other plastomes of Pinaceous genera, the IRs of Pseudolarix and Tsuga are highly reduced, only 449 and 417 bp long, respectively. Their size, gene number, LSC and SSC lengths, and AT content are comparable to those in other Abietoideae genera (table 1). Tsuga and Pseudolarix plastomes share a similar gene repertoire of 35 tRNA genes, four rRNA genes, and 73–74 protein-coding genes (table 1). The Tsuga plastome has one less protein-coding gene than Pseudolarix because its psbI gene is truncated in the type 1 repeat (T1R; see Wu, Lin, et al. 2011) (table 1; supplementary fig. S1, Supplementary Material online). Variations for the total number of genes among Pinaceous genera are due to loss/gain of genes within the T1R. No functional ndh gene has been found in the plastomes of Pseudolarix and Tsuga (supplementary fig. S1, Supplementary Material online), which confirms the loss of all 11 plastid ndh genes from Pinaceae (Braukmann et al. 2009).

Table 1

Comparisons of Plastome Features Among the Five Genera of Abietoideae

Features	Abies koreana	Cedrus deodara	Keteleeria davidiana	Pseudolarix amabilis	Tsuga chinensis
Size (bp)	121,373	119,299	117,720	122,234	120,859
LSC length	66,648	65,052	64,648	65,892	65,105
SSC length	54,197	53,775	52,067	55,444	54,920
Residual IR length	264	426	262	449	417
Pinaceae-specific repeatsa
Type 1 Repeat (T1R)	1,186	1,335	1,286	1,314	1,098
% AT content	61.8	60.9	61.4	61.5	61.9
Total number of genes	113	114	113	113	112
Number of protein-coding genesb	74	75	75	74	73
Number of tRNA genes	35	35	34	35	35
Number of rRNA genes	4	4	4	4	4
Number of duplicated genes	4	5	5	4	3
Within IR/T1R	1/3	1/4	1/4	1/3	1/2

Pinaceae-specific repeats identified by Wu, Lin, et al. (2011).

All ndh genes have been lost from the plastomes of all Pinaceous genera.

Both Pseudolarix and Tsuga plastomes have the A form gene order and contain a pair of T1Rs. However, the T1Rs differ between Pseudolarix and Tsuga, with the former containing a full length psbI gene and being 1,314 bp in length, which is remarkably longer than that of the latter (1,098 bp). Because repeats longer than 200 bp are effective substrates for homologous recombination (Day and Madesis 2007), the T1Rs of both Tsuga and Pseudolarix might also be capable of triggering homologous recombination. No type 2 or type 3 repeat was detected in the plastomes of both species.

Expansion of the AccD Reading Frame in Tsuga

The accD of Tsuga is 1,257 bp, which is longer than the average accD for other Pinaceous genera (969 ± 5 bp). The sequences of accD are highly conserved among the five representative Pinaceous genera (>75% similarity; see fig. 1), with the exception of an about 300-bp insertion that is unique to Tsuga. This insertion is characterized by 23 repeats of the PD/H amino acids (fig. 1).

. 1.—

Comparisons of accD between Tsuga and other Pinaceous genera. (A) mVISTA similarity plots of accD. Blank areas between 300 and 650 bp indicate an insertion specific to Tsuga. (B) Alignment of amino acid sequences showing PD/H tandem repeats in the insertion specific to Tsuga. Repeats are denoted with blue colored-boxes. The histogram below the aligned sequences indicates the level of sequence similarity.

In conifers, expansion of accD was previously discovered in Taiwania and Cephalotaxus with specific tandem repeats characterized by KKD(EY)CDNNE and SDIEED amino acids, respectively (Yi et al. 2013). Including the PD/H tandem repeats in Tsuga, tandem repeats within accD are diverse among conifers. These repeats might have high turnover rates, resembling those within ycf4 of legumes (Magee et al. 2010). Furthermore, Gurdon and Maliga (2014) suggested that in Medicago truncatula, the tandem repeats within accD are recombinationally active and variable among ten ecotypes. Therefore, the repeat is a good population genetic marker. The repeats we discovered in the accD of Tsuga may also be useful in population genetic study of the genus.

AccD codes for the carboxyl transferase β-subunit of the acetyl-CoA carboxylase protein, which is required in fatty acid synthesis (Sasaki and Nagano 2004) and plays a role in leaf development in tobacco (Kode et al. 2005). Positive selection for accD in some angiosperms was proposed to be associated with adaption to various environments (Hu et al. 2015). In Tsuga, we also detected positive selection (dN/dS = 79.04549, P = 0.0319) in the 3’ region of accD where the catalytic sites are located (supplementary fig. S2, Supplementary Material online; Lee et al. 2004); however, its impacts on the evolution of Tsuga require further evaluation.

Phylogeny of Ten Pinaceous Genera Revisited

Overall, 21 taxa were used in the phylogenetic analyses (table 2). Both maximum likelihood (ML) and Bayesian inference (BI) trees have an identical topology (fig. 2), with almost all nodes being strongly supported with 100% bootstrap supports (BS) and 1.0 posterior possibility (PP), except for the trichotomy among Pinus, Picea, and Cathaya. The placement of Cathaya has been inconsistent among many studies; some placed it as sister to Picea (Wang et al. 2000; Lu et al. 2014) and others as sister to Pinus (e.g., Lin et al. 2010). Hence, it is best to regard the three closely related genera as a trichotomy (Nkolongo and Mehes-Smith 2012). Incorporating additional genes from either nuclear or mitochondrial genomes may resolve the trichotomy.

Table 2

Plant Materials and GenBank Accession Numbers of the Plastid Genome Sequences used in the Phylogenetic Analysis

Species of interest		Collection locality	GenBank accession no.	Voucher informationa
Pinaceae
	Pinus koraiensis Siebold & Zucc.	–	AY228468	–
	Pinus thunbergii Parl.	–	NC_001631	–
	Cathaya argyrophylla Chun & Kuang	Sanzhi District, Taiwan	AB547400	Chaw 1486 (HAST)
	Picea morrisonicola Hayata	Xitou Nature Education Area, Taiwan	AB480556	Chaw 1484 (HAST)
	Larix decidua Mill.	Yangmingshan National Park, Taiwan	AB501189	Chaw 1485 (HAST)
	Pseudotsuga sinensis var. wilsoniana (Hayata) L. K. Fu & Nan Li	Wuling Farm, Taiwan	AB601120	Chaw 1487 (HAST)
	Abies koreana E. H. Wilson	Jeju Island, South Korea	KP742350	KHB1465044 (KH)
	Keteleeria davidiana (Bertrand) Beissner	Academia Sinica, Taiwan	AP010820	Chaw 1482 (HAST)
	Tsuga chinensis (Franch.) Pritzel ex Diels.	Taipingshan Forest Park, Taiwan	LC095866	Chaw 1494 (HAST)
	Pseudolarix amabilis (J.Nelson) Rehder	Sanzhi District, Taiwan	LC095867	Chaw 1495 (HAST)
	Cedrus deodara (Roxb.) G.Don	Xitou Nature Education Area, Taiwan	AB480043	Chaw 1483 (HAST)
Araucariaceae
	Araucaria heterophylla (Salisb.) Franco	University of Adelaide, Australia	KM067155	EB1024 (ADU)
	Agathis dammara (Lamb.) Rich.	National Taiwan University, Taiwan	AB830884	Chaw 1490 (HAST)
Podocarpaceae
	Nageia nagi Thunb. O. Kuntze	Academia Sinica, Taiwan	AB830885	Chaw 1491 (HAST)
	Podocarpus totara G.Benn. ex D.Don	–	KC306742	–
Taxaceae
	Amentotaxus formosana H.L. Li	Dr. Cecilia Koo Botanic Conservation Center, Taiwan	AP014574	Chaw 1493 (HAST)
	Cephalotaxus wilsoniana Hayata	Xitou Nature Education Area, Taiwan	AP012265	Chaw 1492 (HAST)
Cupressaceae s.l.
	Cunninghamia lanceolata (Lamb.) Hooker	Longshan Forest Farm, China	KC427270	–
	Juniperus scopulorum Sarg.	–	KF866299	Adams 13594 (BAYLU)
Ginkgoaceae
	Ginkgo biloba L.	Academia Sinica, Taiwan	AB684440	Chaw 1488 (HAST)
Cycadaceae
	Cycas taitungensis Shen, Hill, Tsou & Chen	National Taiwan University, Taiwan	AP009339	Chaw 1489 (HAST)

KH = Korea National Arboretum, South Korea; ADU = The University of Adelaide, Australia; HAST = Herbarium, Biodiversity Research Center, Academia Sinica, Taipei, Taiwan; BAYLU = Baylor University, United States.

. 2.—

Phylogenomic trees inferred from 65 plastid protein-coding genes. The tree framework is based on the ML tree. Values at nodes are arranged by bootstrap supports (BS) for RAxML and posterior probabilities (PP) for MrBayes. Only nodes with support less than 100% BS/1.0 PP are shown. Inset showed the high support of Pinus-Cathaya sisterhood (0.95 PP) based on the BI phylogeny.

Pseudolarix and Tsuga exhibited similar branch lengths to other Pinaceous genera (fig. 2), which generally have slower substitution rates than cupressophytes (Wu and Chaw 2015). Recent molecular studies (e.g., Lin et al. 2010; Lockwood et al. 2013; Lu et al. 2014) and the present results (fig. 2) congruently suggest two separate groups in Pinaceae; one is Abietoideae comprising Abies, Cedrus, Keteleeria, Pseudolarix, and Tsuga; the other consists of all non-Abietoideae genera, including Pinus, Cathaya, Picea, Pseudotsuga, and Larix. This molecular division agrees with the morphological studies of Van Tieghem (1891) and Price et al. (1987), who divided Pinaceae into Abietoid (Cédrées) and Pinoid (Pinées) groups.

Pseudolarix and Tsuga Are Sisters

In figure 2, Cedrus is the only genus that sisters to the clade of the other four Abietoideae genera (i.e., Abies, Keteleeria, Pseudolarix, and Tsuga). The sisterhood of Cedrus and other Abietoideae genera is fully supported (100% BS and 1.0 PP in fig. 2). The presumed alternative relationships, including Cedrus as sister to the other nine Pinaceous genera, the Abies—Keteleeria subclade, and the Tsuga—Pseudolarix subclade, were all rejected by the AU tests (supplementary table S1, Supplementary Material online). Thus, our data reaffirm the position of Cedrus as sister to the remaining Abietoideae genera (Gernandt et al. 2008; Lin et al. 2010; Lu et al. 2014), rather than to the other Pinaceous genera (Wang et al. 2000).

Our two trees congruently indicate the divergence of the other four genera of Abietoideae into two subclades: (1) Abies—Keteleeria and (2) Pseudolarix—Tsuga. Close sisterhood relationships between and within the two subclades received the maximal support (100% BS and 1.0 PP in fig. 2). The likelihood of alternative relationships previously proposed by other studies, such as the (Tsuga,(Pseudolarix,Keteleeria)) suggested by morphological studies (Frankis 1988; Farjon 1990) and (Tsuga,(Pseudolarix,(Abies,Keteleeria))) inferred from cladistics analysis (Hart 1987) and single-copy nuclear genes (Lu et al. 2014), were statistically different on the AU test (supplementary table S1, Supplementary Material online). Therefore, our work clearly supports the sisterhood of Pseudolarix and Tsuga and disagrees with any other alternatives.

Loss of PsaM as A Synapomorphy of Pseudolarix and Tsuga

Figure 3 shows a comparison of residual IRs among the Abietoideae genera. We reannotated the residual IR of Cedrus from 236 (Lin et al. 2010) to 426 bp. Excluding Abies and Keteleeria, the remaining three genera of Abietoideae have residual IRs that include 3’psbA and trnI-CAU. This suggests that the common ancestor of Abies and Keteleeria has shortened its residual IRs to exclude 3’psbA.

. 3.—

Characteristics of the residual IR and T1R in Abietoideae genera. (A) Comparisons of the residual IRs showing the IR contraction to exclude 3’psbA in Abies and Keteleeria. The residual IRs are denoted by grey bars. (B) Specific loss of psaM from the T1Rs of Pseudolarix and Tsuga. The T1Rs of each species are depicted by the green arrows with their lengths indicated. The presence or absence of psaM is marked with blue solid or dashed lines, respectively. Asterisks indicate intron-containing genes. (C) Simplified phylogeny depicting support of a sisterhood relationship between Pseudolarix-Tsuga and Abies-Keteleeria based on characteristics specified in (A) and (B).

In the Abietoideae genera, the T1Rs vary from 1,098 to 1,335 bp (table 1; fig. 3). The T1R in Tsuga is the shortest, containing only partial psbI and lacking psaM. In contrast, although the T1R of Pseudolarix is the second longest, it also lacks psaM. Apparently, loss of psaM is a synapomorphic character inherited from the common ancestor of Tsuga and Pseudolarix (fig. 3). Collectively, these data indicate that the characteristics of the residual IR and T1R have evolved phylogenetically, rather than randomly.

Conclusions

We reaffirm that the three common plastomic characters, i.e., short residual IRs, presence of a T1R, and loss of all plastid ndh genes, signify the plastomes of all Pinaceous genera. The A form observed in both sequenced plastomes also reinforces Wu, Lin, et al. preposition (2011) that the A form is the most primitive among Pinaceae plastomes. We discovered a unique insertion of PD/H tandem repeats that resulted in the expansion of accD in Tsuga, despite the underlying cause remains unclear. In addition, our plastid phylogenomics supports that the five genera of Abietoideae are monophyletic and that they split into (1) Cedrus alone and (2) two sister subclades, Pseudolarix—Tsuga and Abies—Keteleeria.

Material and Methods

Young leaves from T. chinensis (voucher Chaw 1494) and Pseudolarix amabilis (voucher Chaw 1495) were collected from Taipingshan Forest Park and Sanzhi District, Taiwan, respectively (table 2). Voucher specimens were deposited in the herbarium of Biodiversity Research Center, Academia Sinica, Taipei (HAST). Total DNA was extracted following the CTAB protocol (Stewart and Via 1993). P. amabilis was subjected to long-range polymerase chain reaction (PCR) following the protocol in Lin et al. (2010). T. chinensis was sequenced at Yourgene Bioscience (New Taipei City) using the Illumina GAII platform, producing 1 Gb of 100-bp paired-end reads.

Raw reads from T. chinensis were trimmed and de novo-assembled by using the CLC Genomics Workbench v5.5.1 (CLC Bio, Aarhus, Denmark). Contigs < 1 kb and <50× coverage were discarded. Plastome contigs were searched by using the blastn against the K. davidiana plastome with a threshold of E-value < 10−10. Gaps between plastome contigs were closed with PCR using specific primers. The plastome of P. amabilis was assembled from 12 partially overlapping amplicons of 8–16 kb with >8× coverage. The complete plastome sequences were then annotated by using DOGMA (Wyman et al. 2004) and tRNAscan-SE 1.21 (Schattner et al. 2005) with default options. Plastome maps were drawn by using OGDraw (Lohse et al. 2013).

For phylogenetic analyses, a total of 21 taxa were sampled. The collection sites and GenBank accession information are provided in table 2. All protein-coding genes were extracted from the plastomes of 21 taxa and aligned using MUSCLE (Edgar 2004) implemented in MEGA6 (Tamura et al. 2013) with the Align Codons option and default parameters. We used SequenceMatrix (Vaidya et al. 2011) to concatenate the 65 protein-coding genes (supplementary table S2, Supplementary Material online) common to Pinaceae and selected outgroups. The ML and BI trees were constructed from the concatenated matrix by using raxmlGUI v1.3.1 (Silvestro and Michalak 2012) and MrBayes (Huelsenbeck and Ronquist 2001), respectively. ML analysis was conducted with the GTRGAMMA model, which was recommended by jModelTest v2.1.7 (Darriba et al. 2012). The node supports in the ML tree were estimated with 1,000 bootstrap replicates under a majority-rule consensus. The BI tree was evaluated under the GTRGAMMAI model suggested by MrModeltest v2.3 (Nylander 2004). The analysis was run for 1,000,000 generations and sampled every 100 generations, yielding 10,000 trees. The first 25% of trees were discarded as burn-in, and the remaining trees were used to estimate the Bayesian posterior probabilities.

The Approximately Unbiased (AU) test implemented in CONSEL 0.2 (Shimodaira and Hasegawa 2001) was used to assess the probability of alternative relationships among some discordant nodes. We used mVISTA (Frazer et al. 2004) to compare the sequence variability of accD among the representative Pinaceous genera with Tsuga as the reference. Tandem repeats were manually identified in regions of low similarity in the mVISTA plot. Positive selection of accD was detected by using CODEML of pamlX (Xu and Yang 2013) with the branch-site model A (Yang 2007). The branch leading to Tsuga was specified as a foreground branch and the likelihood values of the alternative and null models were calculated by using the options of seqtype = 1, runmode = 0, CodonFreq = 2, model = 2, NSsites = 2, omega = 1 and either fix_omega = 0 (alternative model) or fix_omega = 1 (null model). The likelihood ratio test (LRT) was used to test model fit.