{'fraction_non_alphanumeric': 0.0350709345871033, 'fraction_numerical': 0.02178607549954897, 'mean_word_length': 5.17972972972973, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thirteen New Plastid Genomes from Mixotrophic and Autotrophic Species Provide Insights into Heterotrophy Evolution in Neottieae Orchids

Abstract Mixotrophic species use both organic and mineral carbon sources. Some mixotrophic plants combine photosynthesis and a nutrition called mycoheterotrophy, where carbon is obtained from fungi forming mycorrhizal symbiosis with their roots. These species can lose photosynthetic abilities and evolve full mycoheterotrophy. Besides morphological changes, the latter transition is associated with a deep alteration of the plastid genome. Photosynthesis-related genes are lost first, followed by housekeeping genes, eventually resulting in a highly reduced genome. Whether relaxation of selective constraints already occurs for the plastid genome of mixotrophic species, which remain photosynthetic, is unclear. This is partly due to the difficulty of comparing plastid genomes of autotrophic, mixotrophic, and mycoheterotrophic species in a narrow phylogenetic framework. We address this question in the orchid tribe Neottieae, where this large assortment of nutrition types occurs. We sequenced 13 new plastid genomes, including 9 mixotrophic species and covering all 6 Neottieae genera. We investigated selective pressure on plastid genes in each nutrition type and conducted a phylogenetic inference of the group. Surprisingly, photosynthesis-related genes did not experience selection relaxation in mixotrophic species compared with autotrophic relatives. Conversely, we observed evidence for selection intensification for some plastid genes. Photosynthesis is thus still under purifying selection, maybe because of its role in fruit formation and thus reproductive success. Phylogenetic analysis resolved most relationships, but short branches at the base of the tree suggest an evolutionary radiation at the beginning of Neottieae history, which, we hypothesize, may be linked to mixotrophy emergence.

Introduction symbiosis and usually exchange their photosynthates for mineral nutrients (Smith and Read 2008). Repeatedly, plants evolved the ability to gain organic compounds also from their mycorrhizal fungi, a mixotrophic nutrition also called partial mycoheterotrophy (Selosse and Roy 2009;Hynson et al. 2013;T e sitel et al. 2018). We hereafter use the word mixotrophy for these plants that combine photosynthetic and mycorrhizal carbon sources. In many mixotrophic lineages from different families, species eventually lost photosynthesis and became fully mycoheterotrophic (Leake 1994;Merckx, Freudenstein, et al. 2013). Traits of mycoheterotrophy such as chlorophyll loss, reduction in leaf surface area, development of short clumpy roots, or production of reserve-less "dust" seeds have been first studied (Merckx, Mennes, et al. 2013). An increasing amount of molecular data now clarifies the genetic traits linked with the evolution to mycoheterotrophy.
The plastid genomes (or plastomes) of 24 fully mycoheterotrophic species have been sequenced to date (Graham et al. 2017;Petersen et al. 2018). Their comparative analysis has resulted in a model of sequential plastid genome degradation associated with mycoheterotrophy, which parallels what has been observed for heterotrophic parasitic plants (Funk et al. 2007; Barrett and Davis 2012;Barrett et al. 2014;Wicke et al. 2016;Graham et al. 2017). Genes encoding subunits of the nicotinamide adenine dinucleotide H dehydrogenase (NDH)like complex, which regulates excessive electron flow in the plastidial electron transfer chain, are the first to be lost, followed by most photosynthesis-related genes. Rubisco large subunit (rbcL), adenosine triphosphate (ATP)-synthase, and plastid-encoded RNA polymerase (PEP) genes are often retained in the early stage of photosynthesis loss. Eventually, degradation encompasses housekeeping genes involved in plastid translation and other specific functions. This model summarizes well how the plastid genome evolves in mycoheterotrophic and parasitic plant as they become more and more dependent on organic carbon uptake (Wicke et al. 2016; reviewed in Graham et al. [2017] and Wicke and Naumann [2018]).
Due to the few data available, still little is known about what occurs in the evolutionary steps prior to mycoheterotrophy, that is, in mixotrophic species that retain photosynthetic abilities. Whether these species already experience relaxed selective pressure on photosynthesis remains unclear in the general case. This is a challenging question because of the scarcity of lineages that have retained together autotrophic, mixotrophic, and mycoheterotrophic species. Known examples include the genus Burmannia L. (Burmanniaceae; Bolin et al. 2017), the tribe Pyroleae (Ericaceae; Lallemand et al. 2016), and the tribe Neottieae (Orchidaceae; Selosse and Roy 2009;Gonneau et al. 2014). In the leafless genus Corallorhiza Gagnebin (Orchidaceae), where no full autotroph occurs, comparison of plastomes from mixotrophic (green) and mycoheterotrophic (nongreen) species showed very limited gene loss in the former and relaxed selective constraints on photosynthesis and ATP-synthase genes in the latter (Barrett et al. 2014). Consistently, the selective regime does not differ between mixotrophic Corallorhiza and autotrophic outgroups (Barrett et al. 2014). The absence of autotrophic Corallorhiza as reference yet prevents firm conclusions regarding selective regime in mixotrophs. In the orchid genus Cymbidium Sw., C. macrorhizon is leafless but still holds A gene was presumed pseudogenized if a frameshift indel resulted in a reduction in protein length higher than 25% of its original length.
chlorophyll in the stem and fruits and the importance of its photosynthesis for fruiting has recently been demonstrated (Suetsugu et al. 2018). The plastome of this mixotrophic species shows no gene loss except for a few NDH genes and is mostly under purifying selection (Kim et al. 2018). The two genera Corallorhiza and Cymbidium where photosynthesis already undergoes relaxed selection in mixotrophic species (Barrett et al. 2014) call for extensive analyses in other lineages.
We sequenced plastomes of mixotrophs in Neottieae, an orchid tribe mostly found in northern hemisphere forests. Neottieae include 6 genera and about 200 species (WCSP 2018) with various nutrition types: Palmorchis Barb. Rodr. (all putative autotrophs), Neottia Guett. (many mycoheterotrophs and some probable mixotrophs; T e sitelov a et al. 2012; Yagame et al. 2016;Schiebold et al. 2018), Epipactis Zinn (numerous mixotrophs and some autotrophs; Lallemand et al. 2018), Cephalanthera Rich. (some mycoheterotrophs and many mixotrophs), Limodorum Boehm. (all mixotrophs; Girlanda et al. 2006), and Aphyllorchis Blume (all mycoheterotrophs; Roy et al. 2009). Early studies on mycoheterotrophic and mixotrophic Neottieae contributed to our understanding of the ecophysiology and evolution of plants feeding on fungal carbon (Selosse et al. 2002;Gebauer and Meyer 2003;Julou et al. 2005). The transition to mixotrophy in this group is consistently linked to a shift in type of mycorrhizal fungi (supplementary table S1, Supplementary Material online). Autotrophic orchids associate with a polyphyletic group of saprotrophic and endophytic fungi, the rhizoctonias (Dearnaley et al. 2012), whereas mixotrophic Neottieae use ectomycorrhizal fungi from different genera, which simultaneously form mycorrhizas with surrounding trees (Selosse and Martos 2014). For the purpose of this study, we had to assign a given nutrition type to each species investigated. Even if there is a continuum between autotrophy and mixotrophy (Jacquemyn et al. 2017), we decided to consider a species as mixotrophic when it dominantly associated with ectomycorrhizal fungi and showed an enrichment in 13 C compared with surrounding autotrophs, a usual marker of fungal carbon gain (Hynson et al. 2013; see supplementary table S1, Supplementary Material online, for references regarding species nutrition).
Plastome sequencing of some Neottieae species by Feng et al. (2016) has yielded important results: 1) new hypotheses regarding the controversial phylogenetic relationships between Neottieae genera (Bateman et al. 2005;Roy et al. 2009;G orniak et al. 2010;Xiang et al. 2012), 2) patterns of gene loss in mycoheterotrophs matching the model described above, and 3) no relationship between the chlorophyll level (presence or absence based on visual estimation) and the plastome degradation in the leafless Neottia clade. Their sampling, however, included only four out of six Neottieae genera and one mixotrophic species (Cephalanthera longifolia (L.) Fritsch), limiting the phylogenetic resolution and the analysis of plastome evolution during the emergence of mixotrophy and mycoheterotrophy. Here, we sequenced 13 additional Neottieae plastomes, including representatives of the 2 remaining genera, as well as 9 mixotrophic species (table 1). We investigated how gene content and selective pressure are related to nutrition type and carried out a tentative plastomebased phylogenetic reconstruction of Neottieae.
Plastome Structure, Gene Content, and Selective Regime None of the plastomes of the 13 photosynthetic species sequenced displayed loss of large chromosome regions (>5 kb) or changes in synteny compared with related autotrophic Neottieae (Feng et al. 2016). Their sizes and GC contents were similar, with the exception of Neottia cordata and Limodorum abortivum which have smaller genomes as a consequence of NDH gene loss (table 1). A given gene was presumed pseudogenized if a frameshift indel resulted in a reduction in protein length higher than 25% of its original length. It was considered physically lost when it failed to pass the 77% similarity threshold (compared with closest available reference species) used for gene annotation and was thus undetected. Gene loss (pseudogenization or undetected) was restricted to the NDH complex in a few species (table 1): N. cordata (10 genes lost, ndhI with 148 amino acids vs. 182 in autotrophic Neottia species may also be nonfunctional), L. abortivum (all genes) and Epipactis microphylla (4 genes). An exception was cemA (chloroplast envelope membrane protein), which is lost in L. abortivum similarly to some mixotrophic Corallorhiza (Barrett et al. 2014), which may indicate a nonessential photosynthetic function. A 9-amino-acid deletion and a 28-amino-acid insertion occurred respectively in rpoC1 for Cephalanthera rubra and accD for N. cordata. However, we did not find any stop codons and/or frameshift indels in their sequence, and accD is known to be highly variable in length (Kim and Lee 2004;Gurdon and Maliga 2014;Wicke and Naumann 2018), suggesting that both proteins may still be functional in these species.
For PAML analysis, four nested models allowing the calculation of different x values for some subsets of branches were defined: M0, the null model with one x for the whole tree; M1, which allowed x to differ between outgroups (x 0 ) and Neottieae (x N ); M2, same as M1 but allowing a different x for Palmorchis (x p ) and the remaining Neottieae (x N* ; below noted as Neottieae* ¼ Neottieae without Palmorchis); M3, same as M2 but allowing x to differ between autotrophic (x Na ) and mixotrophic Neottieae* (x Nx , see fig. 1B for species considered as mixotrophic). Models were compared by likelihood ratio tests following this scheme: M1 versus M0, M2 versus M1, and M3 versus M2. For all subsets, successive comparisons of the four models of x variation among branches did not show a better fit for one model against another one, except for the photosynthesis genes for M2 versus M1 (table 2). M2 allowed different x for outgroups (x 0 ¼ 0.12), Palmorchis pabstii (x p ¼ 0.15), and the remaining Neottieae (x N* ¼ 0.10; table 2). These results were not supported by the RELAX analysis, which indicated instead intensification of selective constraints for PEP genes comparing Neottieae* and Palmorchis (k ¼ 1.98; P ¼ 0.04) and for housekeeping genes comparing mixotrophic and autotrophic Neottieae* (k ¼ 1.18; P ¼ 0.02; table 2). But one should consider that these results are not adjusted for multiple comparisons.
The same analysis carried out for individual photosynthesisrelated genes showed a variation in selective regime only for rbcL (table 2 for rbcL, other genes shown in supplementary table S2, Supplementary Material online). rbcL data fitted models M2 much better than M1. In M2, x was lower for Neottieae* (x N* ¼ 0.04) compared with outgroups (x 0 ¼ 0.20) and P. pabstii (x p ¼ 0.52; table 2). This was supported by the RELAX analysis, which showed an intensification of selective constraints in Neottieae* compared with Palmorchis on the one hand (k ¼ 5.24; P ¼ 10 À4 ), and in mixotrophic compared with autotrophic Neottieae* on the other hand (k ¼ 10.97; P ¼ 0.04). Thus, rbcL experiences stronger purifying selection in Neottieae*, and more particularly in mixotrophic species. For individual nonphotosynthesis genes, PAML analysis showed that both matK and ycf1 þ ycf2 fitted M3 better than M2, with higher x for autotrophic (x Na ¼ 0.45, resp. 0.75) compared with mixotrophic Neottieae* (x Nx ¼ 0.23, resp. 0.55; table 2). This was, however, not supported by the RELAX analysis (table 2).
To summarize, we did not find evidence for relaxed selection in plastomes of mixotrophic species. By contrast, both PAML and RELAX analyses showed stronger selective pressure in mixotrophic Neottieae* compared with Palmorchis for rbcL.

Plastome Evolution in Mixotrophic Neottieae
The observed loss of NDH genes in three plastomes confirms that this gene complex is sensitive to evolutionary decay in Neottieae (Feng et al. 2016) and more generally in orchids (Barrett et al. 2014;Luo et al. 2014;Kim et al. 2015;Lin et al. 2017;Kim et al. 2018) and some other photosynthetic lineages (Wicke and Naumann 2018). Interestingly, it does not correlate with species nutrition. For example, the 11 genes were lost in L. abortivum and N. cordata (if ndhI is indeed nonfunctional). The first is mixotrophic and highly dependent on fungal carbon because its photosynthesis does not even compensate for its respiration (Girlanda et al. 2006); the second displays a more versatile use of fungal resources (Schiebold et al. 2018), which may even be optional (T e sitelov a et al. 2015). Similarly, E. microphylla is the only species of the mixotrophic Epipactis section (E. albensis, E. atrorubens, E. helleborine, and E. purpurata) (Jin et al. 2014) for which NDH gene loss occurred, whereas some losses have been reported in E. veratrifolia and E. mairei (Feng et al. 2016) that are supposed autotrophic (pending more investigations). This suggests a species-specific decay of NDH genes in orchids, independently of their autotrophic or mixotrophic habit. In addition, we did not observe any evidence for intensification or relaxation of selective pressure on these genes (table 2). We suppose that relaxation on NDH genes occurs in the whole orchid clade, making any difference among orchid species hardly detectable. Changes in nutrition type do not seem to alter this pattern.
Photosynthesis-related genes neither are lost nor experience relaxed selective constraints in mixotrophic species but, in Neottieae*, they showed a slightly lower x; rbcL was even under stronger purifying selection compared with P. pabstii (table 2). Some parts of the plastomes unrelated to photosynthesis also experienced higher selective constraints in mixotrophs, as shown by 1) lower x for matK and ycf1 þ ycf2 in mixotrophic than in autotrophic Neottieae* and 2) the RELAX analyses of PEP and housekeeping genes (displaying k > 1). These results are concordant with what has been reported for mixotrophic Corallorhiza, which 1) retain all photosynthesisrelated genes except cemA and psbM, 2) display apparently unimpaired photosynthesis, and 3) show no difference in plastome selective regime compared with autotrophic references (Barrett et al. 2014). This also agrees with the retention of purifying selection on plastidial genes observed in the mixotrophic C. macrorhizon (Kim et al. 2018). However, our results contrast with other kinds of mixotrophy where relaxed purifying selection has been observed for photosynthesis, PEP and ATP-synthase genes, namely 1) in carnivorous plants (psa, atp, and rpo genes in Wicke et al. [2014]) and 2) in obligate hemiparasites (psa, psb, pet, atp, and rbcL in Petersen et al. [2015]; photosynthesis and atp genes in Wicke et al. [2016]). For this latter trophic strategy, retention of selective constraints on photosynthesis genes has nevertheless been reported for photosynthesis genes in Cuscuta (McNeal et al. 2007) and also for some housekeeping genes in Viscum (accD, cemA, clpP, and ycf2, Petersen et al. 2015). However, these comparisons between studies have to be taken with caution given the differences in scales among them or in a given study (individual genes vs. gene subsets).
In Neottieae, depending on the plastid genes considered, mixotrophic species thus feature either no change or a lowly supported increase in selective pressure compared with other mixotrophic groups such as Corallorhiza, carnivorous or parasitic plants ( Barrett et al. 2014;Wicke et al. 2014Wicke et al. , 2016Petersen et al. 2015). Photosynthesis therefore seems to remain advantageous in mixotrophic Neottieae, and this is congruent with physiological evidence for the importance of photosynthetic carbon in fruit development and reproduction in these species (Roy et al. 2013;Bellino et al. 2014;Gonneau et al. 2014;Lallemand et al. 2019). We hypothesize that the retention of genes and purifying selection in the plastome of C. macrorhizon (Kim et al. 2018) indicates a similar crucial role of photosynthesis in reproduction, as illustrated by Suetsugu et al. (2018).

Neottieae Phylogeny and Mycoheterotrophy Evolution
Phylogenetic analysis of the 28 available plastomes confirmed Neottieae monophyly with Palmorchis in sister position to the remaining species ( fig. 1A). Neottia, Epipactis, and Cephalanthera genera were each monophyletic and most intrageneric relationships were well defined. Aphyllorchis and Limodorum were sister genera but with mild support from maximum likelihood (ML) analysis (bootstrap ¼ 74, fig. 1A). However, very short branches with low support prevented strong conclusions about intergeneric relationships ( fig. 1A). Because mycoheterotrophic Neottia species display accelerated evolutionary rates (see long branches in fig. 1A), we removed them to increase alignment quality: We then obtained a congruent tree with better support for most but not all clades (supplementary fig. S1A, Supplementary Material online). In particular, Neottia, Epipactis, Aphyllorchis, and Limodorum clustered together in a monophyletic group sister to Cephalanthera with better support from ML analysis (bootstrap ¼ 87 vs. 71). When restricting the analysis to the plastidial coding sequences (CDS) of nonmycoheterotrophic species, that is, excluding noncoding regions to allow better alignment, Limodorum was sister to Epipactis but Neottia no longer clustered with these genera in a monophyletic group (supplementary fig. S1B, Supplementary Material online). Combining these results suggests that at least Cephalanthera is sister clade to the rest of Neottieae. Analysis based on ITS (internal transcribed spacer of nuclear ribosomal DNA, including ITS1-5.8S rRNA-ITS2) confirmed Aphyllorchis and Limodorum as sister genera but failed to clarify further intergeneric relationships (supplementary fig. S1C, Supplementary Material online).
Phylogenetic relationships between the Neottieae genera have long been controversial and studies using different species and marker combinations yielded contradictory results, without any clear intergeneric relationships (ITS and the plastid trnL intron in Bateman et al. [2005]; ITS, trnS-G spacer, and rbcL in Roy et al. [2009]; ITS, Xdh, rbcL, matK, psaB, and trnL-F spacer in Xiang et al. [2012]; ITS, rbcL, and matK in Zhou and Jin [2018]). More confusingly, the phylogeny of Feng et al. (2016), based on a subset of the plastomes we analyzed, disagrees with our results, displaying Cephalanthera and Aphyllorchis as sister genera, and Epipactis and Neottia as successive sisters. This inconsistency may result from the difference in the number and diversity of species included in the analysis. We added 13 species and 2 new genera (Palmorchis and Limodorum) but failed to reproduce their results when considering only their subset of species, suggesting that this is rather due to different phylogenetic methods and/or to evolutionary signals conveyed by the plastome regions used. Feng et al. (2016) selected conserved plastome regions, whereas we made a dual analysis of whole plastome alignment and CDS-restricted alignment. Moreover, we chose not to filter or manually edit alignments following Tan et al. (2015) for accurate phylogenetic inferences (see also Hallas et al. 2017).
Neottieae intergeneric relationships thus remain difficult to estimate confidently. Our plastome-based analysis showed that the branches separating these genera are very short (<10 À3 substitutions per site). Thus, the common ancestor of Neottieae* (i.e., Neottieae with exclusion of Palmorchis) likely experienced rapid speciation, with little time for informative mutations to accumulate (Glor 2010). Such a fast evolution of ancestral lineage makes it prone to hybridization and incomplete lineage sorting (Glor 2010;Wang et al. 2014), so that independent loci (e.g., nuclear and plastidial markers) may encapsulate different phylogenetic stories and make any combined analysis unreliable. Another possibility is heterotachy (Fitch and Markowitz 1970), a variation (here a decrease) in evolutionary rate in the ancestral lineage of Neottieae*. Arbitrating between fast speciation and heterotachy would be particularly challenging given the very short branch lengths. Confident reconstruction of early speciation in Neottieae will be complicated, but different gene trees could provide interesting insights into early Neottieae evolution and possible hybridization events.
We also wanted to compare two possible evolution of nutrition in Neottieae*. Based on our plastid phylogeny, we display two hypotheses, namely 1) an autotrophic ancestor and 2) a mixotrophic ancestor for Neottieae* (see Materials and Methods for reconstruction of intermediate characters). Assuming an autotrophic common ancestor, mixotrophy evolved at least four times ( fig. 1B, gray boxes above the branches). Conversely, assuming a mixotrophic common ancestor would give three reversions to autotrophy ( fig. 1B; white boxes below the branches). Each scenario implies two emergences of mycoheterotrophy, but unfortunately, the lack of information about the trophic status of some taxa (such as E. mairei or Neottia fugongensis) and the absence of several mycoheterotrophic clades limit our analysis. For example, green Neottia other than N. cordata or Epipactis species outside the Epipactis section may actually show mixotrophic nutrition. Assuming a mixotrophic common ancestor implies one change less than an autotrophic one but includes reversion to autotrophy. Although unlikely at first glance, reversion is allowed because mixotrophy is flexible and displays a continuum from autotrophy to mycoheterotrophy (Jacquemyn et al. 2017). Most importantly, the selective pressure maintaining photosynthesis genes reported here allows a reversion from mixotrophy to autotrophy, because potential for autotrophy remains intact in mixotrophic Neottieae. The question whether mixotrophy is a plesiomorphy of Neottieae thus remains open, pending for analyses of more species.
A final speculation comes at this point: If mixotrophy turns out to be indeed the ancestral state of Neottieae, could it have been a key factor in the fast early diversification speculated above? Mixotrophy allows plants to grow in light-limited forest understory and so to occupy a broad new ecological niche (T e sitel et al. 2018). Orchids show important dispersion abilities and most mixotrophic Neottieae do not seem very selective regarding the fungi they associate with (T e sitelov a et al. 2012). One can therefore speculate that an early mixotrophic Neottieae lineage could have rapidly colonized distant and heterogeneous forest patches, with ability to keep or loose photosynthetic abilities, thus favoring radiation.

DNA Extraction and Sequencing
Total DNA was extracted from the leaves or fruits of the following Neottieae species (table 1 and supplementary table S1 Porebski et al. (1997).
Given low yield due to DNA degradation in some old samples, the Accel-NGS 1S Plus DNA Library Kit for Illumina with the 1S Plus Indexing kit (Swift Biosciences Inc., Ann Arbor, MI) was used for preparation of libraries, following the manufacturer's instructions. We did so for the set of samples handled in EPGV laboratory (supplementary table S3 Read Cleaning, Plastome Assembly, and Annotation Redundant reads were removed with a homemade Cþþ script. The redundancies were sought by pairs read1-read2, the redundancy with the best Phred score was kept. Trimming of the low-quality 3 0 end was done using a homemade Cþþ script, one base after the other until reaching a base with a Phred score > 30 or until the mean Phred score of the read was above 30. Pairs with at least one read shorter than 30 bp after this step or with N were removed. Read number before and after cleaning are indicated in supplementary table S4, Supplementary Material online. For samples that were prepared with the Accel-NGS 1S Plus DNA Library Kit, adapter trimming was conducted using BBDuk Trimmer plugin version 37.64 in Geneious version 11.1.3 (http://www.geneious.com; last accessed September 2018; Kearse et al. 2012) with the following parameters: adapters ¼ Illumina Truseq DNA adapters, ktrim ¼ r, k ¼ 27, hdist ¼ 1, edist ¼ 0, mink ¼ 6, minlength ¼ 10, trimbyoverlap ¼ t, and minoverlap ¼ 6. Addition of a low complexity tail to the 3 0 end of fragments was part of the library preparation kit and had to be removed from the beginning of reads 2. Following the manufacturer's recommendations, trimming an additional 10 bp at the 5 0 end was done for all reads.
De novo assembly was done differently depending on the species (supplementary table S3, Supplementary Material online). 1) A subset of 25% of the reads was assembled using the Geneious algorithm with medium-low sensitivity parameters. Resulting contigs were mapped on the plastome of Neottia ovata (L.) Bluff & Fingerh. (KU551271), Cephalantera longifolia (L.) Fritsch (KU551263), or E. veratrifolia Boiss. & Hohen (KU551267), depending on the species. Contigs corresponding to the plastome were then dissolved and reassembled with medium sensitivity parameters to increase assembly quality. The resulting contigs were again mapped to the reference genome and large single copy (LSC), inverted repeat B (IRB), and small single copy (SSC) regions were isolated. IRB was duplicated and reverse complemented to obtain inverted repeat A (IRA). Concatenation of LSC, IRB, SSC, and IRA led to the final plastome sequence.
2) The CLC Genomics Workbench de novo assembly procedure was used. In this latter case, reads were first trimmed with the following settings: quality limit ¼ 0.01 (corresponds to Q ! 20), remove adapters ¼ yes (Illumina Truseq adapters), discard short reads ¼ yes (limit 50 nt for MiSeq data, 25 nt for HiSeq data). Then de novo assembly was performed with the following settings: automatic bubble size ¼ yes, minimum contig length ¼ 1,000 bp, automatic word size ¼ yes, perform scaffolding ¼ yes, and auto-detect paired distances ¼ yes. Contigs of plastid origin were selected based on the results of BLAST search with Phalaenopsis aphrodite plastome used as query. Contigs corresponding to LSC, IR, and SSC were concatenated accordingly and resulting sequences were verified and corrected by performing the gap closure method using back mapping of the reads. IRs were duplicated and reverse complemented to obtain the final parts (IRA) of the plastomes.
Coverage values of the contigs used for assembly are indicated in supplementary table S4, Supplementary Material online. Depending on the species, the phylogenetically closest reference plastome was used for annotations transfer in Geneious with a 77% similarity threshold that turned out, after preliminary attempts, to be low enough to detect similar gene sequences and sufficiently high to avoid false positives. Manual correction and addition of missing information was then carried out for each plastome. Annotated plastome sequences were submitted to GenBank under the accession numbers MH590345-MH590357 (table 1).

Phylogenetic Inferences
Including the 12 already published sequences, we ended up with a set of 25 Neottieae plastomes and added 3 outgroup species from the same subfamily Epidendroideae (supplementary table S1, Supplementary Material online). The four mycoheterotrophic Neottia were either kept or removed to test whether their rapidly evolving plastomes (much higher than the two other mycoheterotrophs; fig. 1A) and the increased length and quality of the alignment obtained without them alter phylogeny reconstruction. The second copy of the inverted repeat region (IRA) was removed from all plastomes and the resulting sequences were aligned a first time using the Mauve progressive algorithm, as implemented in Geneious (parameters set to automatically calculate seed weight and minimum locally collinear block (LCB) score, compute LCBs, full alignment), to detect potential rearrangements (Darling et al. 2004). The LCBs obtained were individually aligned a second time using the MAFFT online service (Katoh et al. 2017) with automatic selection of alignment strategy, gap opening penalty set to 3.0, offset value set to 0.1 and other parameters left as default. The resulting alignments were then concatenated. Model selection was done with the ModelFinder algorithm (Kalyaanamoorthy et al. 2017) implemented in the IQ-TREE program (Nguyen et al. 2015). The model with lowest Bayesian information criterion score was used for ML analysis with IQ-TREE (Nguyen et al. 2015). The first model with the lowest Bayesian information criterion compatible with MrBayes v3.2 parameters (Ronquist et al. 2012) was chosen for Bayesian inference with this program. A complementary analysis using the nuclear marker ITS, including ITS1-5.8S rRNA-ITS2, was done using available Neottieae sequences (supplementary  table  S1, Supplementary Material online).
On the plastid phylogeny ( fig. 1B), we displayed two possible evolution of nutrition based on the hypotheses of 1) an autotrophic ancestor or 2) a mixotrophic ancestor for Neottieae*. Intermediate characters for other nodes were proposed in order to minimize the number of transitions between trophic types (autotrophic, mixotrophic, or mycoheterotrophic), that is, the most parsimonious scenario under each hypothesis, considering the nutritional status of extant species (supplementary table S1, Supplementary Material online). A given species was considered as mixotrophic when it dominantly associated with ectomycorrhizal fungi and showed an enrichment in 13 C compared with surrounding autotrophs, a usual marker of fungal carbon gain (Hynson et al. 2013; see supplementary table S1, Supplementary Material online, for references regarding species nutrition).

Selective Regime Analyses
All CDS were extracted from the same set of plastomes used for phylogeny reconstruction, except the six mycoheterotrophic species. Except for NDH complex, genes that were lost or pseudogenized in at least one of the remaining species were discarded from subsequent analyses (see table 2 for the remaining genes). For NDH, because some species had lost all the genes, we rather only kept for subsequent analysis only the 15 species, which retained the 11 NDH genes (table 2). CDS were aligned based on their amino acid sequences with subsequent backtranslation using the AlignTranslation function (Wright 2015) of the DECIPHER package (Wright 2016) in the R environment for statistical computing (R Development Core Team 2007). Functional gene groups were built by concatenating CDS for photosynthesis, NDH, ATP-synthase, PEP, and housekeeping genes (table 2). ycf1 and ycf2 were combined and analyzed independently because of their specific evolutionary rates (Barrett et al. 2014). In addition, some genes were analyzed individually when long enough for confident parameter estimation: ccsA, psaA, psaB, psbA, psbB, psbC, psbD, rbcL, accD, and matK (see supplementary table S2, Supplementary Material online). x, the ratio of nonsynonymous substitutions per nonsynonymous site (d N ) to synonymous substitutions per synonymous site (d S ), was used to estimate selective pressure for different sets of genes. The codeml program implemented in PAML version 4.9 (Yang 2007) was used to compute x values, with the branch model and the tree shown in figure 1B. The codon frequency parameter was set to F3X4. Four nested models allowing the calculation of different x values for some subsets of branches were defined: M0, the null model with a unique x for the whole tree; M1, which allowed x to differ between outgroups (x 0 ) and Neottieae (x N ); M2, same as M1 but allowing a different x for Palmorchis (x p ) and the remaining Neottieae (x N* ); M3, same as M2 but allowing x to differ between autotrophic (x Na ) and mixotrophic Neottieae* (x Nx , see fig. 1B for species considered as mixotrophic). Models were compared by likelihood ratio tests following this scheme: M1 versus M0, M2 versus M1, and M3 versus M2. The statistic of the test (K) was calculated as two times the difference in model likelihood and compared with a chi-squared distribution with degrees of freedom equaling the difference in parameters number between the two models. P values were obtained using the function pchisq of the R stats package.
We completed the PAML analysis with a specific test dedicated to the detection of changes in selective pressure (relaxation or intensification). This was carried out with the program RELAX, available at the Datamonkey webserver http://datamonkey.org/relax; last accessed March 6, 2019. RELAX allows more subtle estimation of x variations by considering different categories of sequence sites that can evolve under different selective constraints (Wertheim et al. 2015). x distributions are estimated for a test (T) and a reference (R) subset of branches and a selection intensity parameter k such that x T ¼ x k R is introduced. The program then compares the goodness of fit of the data to a branch-site evolutionary model with k ¼ 1 (null model) or k being a free parameter (alternative model). When the alternative model shows a better fit following a likelihood ratio test, k < 1(resp. k > 1) indicates a relaxation (resp. an intensification) of selective constraints in the test group. Three different tests were done with three different designs of test versus reference groups, copying what has been done with the PAML nested analysis (see supplementary fig. S2, Supplementary Material online, for group designs).

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online.

Acknowledgments
We acknowledge Aur elie Chauveau for preparation of the libraries, Elodie Marquand and Aur elie Canaguier for data processing, two anonymous referees and Dennis Lavrov for their corrections, and David Marsh for English corrections. This work was funded by the Fondation de France (funds from the Fondation Ars Cuttoli & Paul Appell, project name : EVOLUTION ET PHYSIOLOGIE DES PLANTES QUI SE NOURRISSENT DE CHAMPIGNONS (MYCOHETEROTROPHIE)) and the National Science Center (Poland; grant no. 2015/18/A/NZ8/00149). We thank the Commissariat a l' Energie Atomique et aux Energies Alternatives -Institut de G enomique/Centre National de G enotypage for conducting the DNA quality control and for giving the INRA-EPGV group access to their Illumina Sequencing Platform. Plastomes of Ce. damasonium and L. abortivum were sequenced with support of budgetary subsidy to IITP RAS (Laboratory of Plant Genomics). We thank the MSU Herbarium (MW) for providing a sample of Ce. longibracteata Blume MW0047886 and "Service de Syst ematique Mol eculaire" (UMS2700 MNHN/CNRS) for granting access to its technical platform.