Leaf tissue was obtained from 17 individuals in total: 11 accessions were obtained from the Atlanta Botanical Garden’s living conservation collection (S. oreophila, S. jonesii, S. alata, S. alabamensis and S. rubra) and six accessions were obtained from two field sites (S. rubra subsp. rubra and S. rubra subsp. viatorum). DNA was extracted from silica dried samples using the Qiagen DNeasy Plant Mini Kit. Library prep was performed using the Kapa Biosystems HyperPlus Kit using iTru adapters (Glenn et al., 2019). Libraries were pooled at equal concentrations and enriched for putative single-copy orthologs enrichment using the Angiosperms353 bait set (Johnson et al., 2019). The enriched pool was sequenced on an Illumina NextSeq 500 at the Georgia Genomics and Bioinformatics Core using a High Output 300 cycle flow cell generating 150bp paired-end reads.

In addition, sequencing reads from Stephens et al. (2015) were downloaded from NCBI Short Read Archive. The Stephens et al. data set includes 71 accessions of Sarracenia and 4 accessions of outgroups in Sarraceniaceae (Heliamphora minor and Darlingtonia californica).

All raw reads were trimmed using Trimmomatic (v. 0.39) (Bolger et al., 2014). Both the new data set and the data set obtained from Stephens et al. were sequenced from libraries enriched for targeted nuclear loci. However, the majority of the reads from both data sets are off-target. Stephens et al. reported an average of 1.6% of reads on target, and analysis of the new data set revealed that less than 1% of the reads were on target. The large proportion of off-target reads enable the assembly of the plastome.

Initial de-novo plastome assembly was attempted with GetOrganelle (v. 1.7.5.2) (Jin et al., 2020). GetOrganelle often produced two assembly versions differing only in the orientation of the short single copy regions (SSC). SSC orientation was determined by aligning assemblies to the reference plastome (Clethra L. delavayi, Genbank accession NC_041129) using MUMmer (v. 4.0.0) (Kurtz et al., 2004), and only the assemblies with concordant SSC orientation were retained.

GetOrganelle did not generate complete de-novo plastome assemblies from every sample. In these cases, the following reference-based pipeline was used. Reads were aligned to one of the complete Sarracenia plastome assemblies using BWA (v. 0.7.17) (Li et al., 2009). The aligned reads were then extracted and assembled de-novo using SPAdes (Bankevich et al., 2012). Afin (https://github.com/afinit/afin) was used to extend the resulting contigs and fuse any contigs with significant overlap. At this stage, assemblies were either mostly complete (1-3 contigs consisting of the large single copy region (LSC), short single copy regions (SSC), and one IR), or they were more fragmented. The mostly complete assemblies were manually pasted together. The IR boundaries were verified by mapping reads to the assemblies and identifying the coordinate where half of the reads spanned the IR and LSC and the other half spanned the IR and SSC.

Complete plastome assemblies were annotated using PGA (Qu et al., 2019). Fragmented assemblies were aligned to one of the complete, PGA annotated plastomes using the Minimap2 (v. 2.17) (Li, 2018) plugin in Geneious. The “transfer annotation” function was used before generating a consensus sequence.

Coding sequences (CDS) from 79 plastid genes were extracted from the annotated assemblies and aligned with MAFFT (v. 7.470) (Katoh and Standley, 2013). All resulting gene alignments were concatenated. Regions of the concatenated alignment that were poorly aligned or had gaps in 50% or more of the samples were filtered out of the gene alignments using Gblocks (v. 0.91b) (Castresana, 2000). A maximum-likelihood phylogeny was estimated from the concatenated gene alignments using IQ-Tree (v. 2.0.6) (Nguyen et al., 2015). 1000 bootstrap replicates were performed using UFBoot (Minh et al., 2013). The GTR + F + R4 substitution model was used.

To differentiate between incomplete lineage sorting (ILS) and chloroplast capture, a tree simulation approach similar to Folk et al., 2017 (Folk et al., 2016) was used. Plastome trees under ILS were simulated using the dendropy python package (v. 4.5.2) (Sukumaran and Holder, 2010) with the species tree from Stephens et al. (2015) as a guide tree. Since plastomes are effectively haploid and inherited uniparentally, plastomes have one quarter of the effective population size of diploid nuclear loci. Since the guide tree used for these simulations was estimated exclusively using nuclear loci, its branch lengths were scaled by four to account for the effective population size differential between plastomes and nuclear loci. A distribution of tree discordance under the null hypothesis of ILS was generated by calculating a tree distance metric [information-based generalized Robinson-Foulds distance (Smith, 2020)] between 1000 simulated trees and the species tree. Then the distance between the empirical plastome tree from this study and the species tree was calculated and compared to the null distribution. Since the empirical plastome tree has samples that are not in the Stephens et al. (2015) species tree, those tips were dropped from the plastome tree to enable calculating distance.

Fourteen complete, circularized plastomes have been assembled and annotated including the following Sarracenia species: S. jonesii, S. alabamensis, S. oreophila, S. rubra subsp. gulfensis, S. rubra subsp. rubra, and S. rubra subsp. viatorum. Average assembly statistics for the all assemblies are shown in Table 1 . The assembly pipeline for fragmented assemblies recovered an average of 114kbp of plastome sequence, almost tripling the 42kbp recovered in Stephens et al. (2015). The use of different references is one potential factor explaining this difference; this study used a complete Sarracenia plastome (Ericales) as a reference whereas Stephens et al. (2015) used a plastome from Vitis vinifera (Vitales). Eighty protein-coding genes were extracted from assemblies, and sequences were aligned for all samples, and alignments were concatenated for the phylogeny estimation.

All complete Sarracenia plastomes include some pseudogenized plastome-encoded genes. With the exception of ndhB and ndhE, all ndh genes either have been pseudogenized due to premature stop codons or large deletions ( Figure 1 ). Similarly, all samples contain a premature stop codon within the rps12 gene.

Consistent with Stephens et al. (2015), no species were found to exhibit monophyly of their plastomes, and the plastid tree is highly incongruent with the published species tree ( Figure 2 ). Support values across the backbone of the tree are all greater than 70, and most internal nodes are highly supported as well ( Figure 2 ). Branch lengths within Sarracenia are generally very short in comparison to the outgroups. An exception is the split at the base of the Sarracenia clade. This branch splits Sarracenia into two distinct plastid lineages. These main lineages are arbitrarily termed clade A and clade B ( Figure 2 ). Clade B contains all sampled individuals of minor, oreophila, jonesii, and purpurea var. montana, and clade A contains all sampled individuals of alata and purpurea (excluding var. montana). All other species are split across these two main lineages (flava, psittacina, rubra, and leucophylla).

The tree distance metric that was used ranges from 0 (an identical tree) to 1 (the most distal tree). The plastome trees simulated under the pure coalescent model have distances from the species tree ranging from 0.29 to 0.56, while the distance from the empirical plastome tree is 0.73 ( Figure 4 ). A T-test using the distribution of simulated plastome tree distances as the null distribution gives a p-value of >2.2e-16, rejecting the null hypothesis of ILS causing the discordance alone.

Independent pseudogenization or complete loss of ndh genes has been shown in many plant lineages, including holoparasitic, hemiparastitic, and carnivorous plant lineages (Barrett et al., 2014; Lin et al., 2017; Cao et al., 2019; Gruzdev et al., 2019; Nevill et al., 2019). Functional ndh genes are rarely found in non-photosynthetic parasitic plants, and the loss of ndh genes is strongly correlated with the transition to heterotrophy in parasitic plant lineages (Wicke et al., 2016). Since plastid encoded ndh genes are thought to optimize photosynthetic chemistry in fluctuating or stressful environments [reviewed in (Sabater, 2021)], the loss of ndh genes in parasitic lineages that are no longer fully dependent on photosynthesis as a source of carbon is unsurprising. In carnivorous plants, however, evidence for significant heterotrophic uptake of carbon is limited (Rischer et al., 2002), and a transition to full heterotrophy seems unlikely, so this line of reasoning does not explain the independent pseudogenization of functional ndh genes across carnivorous plant lineages. It is possible that the acquisition of organic nitrogen has an interaction with photosynthetic chemistry that relaxes the need for ndh. As Nevill et al. (2019) noted, organic nitrogen acquisition bypasses the need to assimilate nitrate using photosynthetically-derived reductant. Alternatively, the pseudogenization of ndh genes in parasitic plants and carnivorous plants could be due to unrelated mechanisms. The pseudogenization of almost all of the ndh genes across the genus Sarracenia shown here provides further evidence that carnivorous plants do not require these genes. Sequencing of full plastomes from other carnivorous species would reveal if the pseudogenization of ndh occurs early in carnivorous plant evolution.

The plastome phylogeny in this study shows a similarly extreme level of discordance with the species tree as that of the Stephens et al. (2015) plastome phylogeny. That study ascribed the discordance to a combination of chloroplast capture and a lack of informative polymorphisms in the chloroplast sequence. A third source of discordance, ILS, is considered here. A lack of informative polymorphisms is not an issue here, as almost all the plastome coding sequences are used and the resulting phylogeny has high bootstrap values across the spine, suggesting that there is sufficient evidence that major clades within the tree are correct.

To distinguish between the two remaining sources of discordance, plastome phylogenies under ILS were simulated. The simulated phylogenies showed much lower levels of discordance with the species tree than the empirically estimated plastome. To simulate the plastome phylogenies, the branch lengths of the guide tree were multiplied by four due to the assumption that the chloroplast is inherited matrilineally in Sarracenia like most seed plants (Mogensen, 1996). Since this assumption hasn’t been empirically proven and biparental inheritance of the chloroplast is possible, simulations with branch lengths multiplied by two were performed and show similar results ( Supplemental Data ).

There is ample signal of introgression in the plastome, but Stephens et al. (2015) reported no evidence of gene flow in the nuclear data. A search through the nuclear gene trees revealed that none of the trees had a similar topology to the plastome tree, but some trees did exhibit a high degree of discordance with the species tree, possibly due to occasional nuclear gene introgression. Cytonuclear discordance is commonly observed and is attributed to introgression in plant and animal systems (Rieseberg and Soltis, 1991; Berthier et al., 2006; Gernandt et al., 2018), including several instances where there is limited signal for introgression in nuclear data (Winkler et al., 2013; Good et al., 2015; Folk et al., 2016; Rose et al., 2020). However, the mechanism for organellar introgressions without accompanying nuclear loci is poorly understood (Rieseberg and Soltis, 1991; Folk et al., 2018). Sarracenia is a genus where hybridization is common and thus some level of nuclear introgression might be expected. The extreme level of chloroplast capture and lack of signal for nuclear gene flow in Sarracenia illustrates the comparative ease of introgression of organelles over nuclear loci.