We collected G. mupinensis and P. tubiflorus which represent two holoparasitic genera from the provinces Sichuan and Heilongjiang, China, respectively. Voucher specimens were deposited in the herbarium of Fudan University (FUS). Total genomic DNA was extracted from fresh tissues dried with silica-gel following the CTAB extraction method (Doyle and Doyle, 1987).

Specific primers of five genes (Table S1) were used for polymerase chain reaction (PCR). Primers used to amplify the nuclear genes PHYA and PHYB as previously described (Bennett and Mathews, 2006; McNeal et al., 2013) did not work in G. mupinensis and P. tubiflorus; we designed new primers (Table S1) according to the sequences of closely related taxa inferred from phylogeny of plastid (rps2 and matK) data. Each reaction volume of 50 μl contained ~150 ng total DNA, 5 μl of 10X PCR buffer, 6 μl of MgCl₂ (2.5 mmol·l⁻¹), 8 μl of dNTP mixture (2.5 mmol·l⁻¹), 3 μl of each primer (10 μmol·l⁻¹), and 2.5 units of Red Taq DNA polymerase. The running programs for these genes are presented in Table S2. PCR products were visualized on 1% agarose gels and purified using Gel Extraction System B Kit (BioDev-Tech, Beijing, China) according to the manufacturer's instructions. The purified products were sequenced directly with BigDye Terminator 3.1 Cycle Sequencing Kit (Applied Biosystems) and run in ABI PRISM 377XL DNA Autosequencer. All sequences generated in this study were submitted to GenBank (accession numbers: KY706614–KY706632).

The sequences of five genes from G. mupinensis and P. tubiflorus were assembled by the software Seqman II 5.05 (DNAStar, London, UK). Sequence data of other species in Orobanchaceae were obtained from GenBank (Table S3). The alignment was undertaken by Clustal X (Larkin et al., 2007) and adjusted manually. All data sets were partitioned by gene, and gaps were treated as missing data. Modeltest 3.5 (Posada and Crandall, 1998) was executed to select the best fitting evolutionary model for each partition under the Akaike Information Criterion (AIC). Maximum likelihood (ML) and Bayesian inference (BI) analyses were performed on the independent and combined data sets for ITS, PHYA, PHYB, rps2, and matK. ML analyses were run in RAxML 7.0.4 (Stamatakis, 2006); the general time reversible (GTR) model and gamma distribution were used for all partitions (McNeal et al., 2013) because it is impossible to specify different models to different partitions in RAxML. Support values for branches of the ML trees were assessed based on 1,000 bootstrap replicates. Before combining all data sets, the partition homogeneity test was used to check the congruence among the five-gene data sets in PAUP^* 4.0b10 (Swofford, 2002). Bayesian analyses were performed in MrBayes 3.1.2 (Huelsenbeck and Ronquist, 2001). Each partition was assigned its own nucleotide substitution model. We ran 25 million generations sampling every 1,000 generations. Stationarity of each analysis was assessed using Tracer 1.6.0 (Rambaut and Drummond, 2012) by checking the effective sample size (ESS) values. Completion was determined when the average standard deviation of split frequencies fell below 0.01. The 50% majority rule consensus trees were obtained after the first 25% of the samples were removed as burn-in.

In order to compare the resolution of different partitioning schemes, we also divided the data by codon to construct the phylogeny of Orobanchaceae based on single and combined gene data, respectively. Each protein-coding gene (PHYA, PHYB, rps2, matK) was partitioned by codon (1st, 2nd, and 3rd codon position) and the best fitting models were selected by Modeltest 3.5. The concatenated five-gene data set was partitioned into 13 blocks, including one for the ITS region and 12 for the 1st, 2nd, and 3rd codon positions of each protein-coding gene. ML analyses were conducted with RAxML 7.0.4 under the GTR + G model and Bayesian inference was run in MrBayes 3.1.2 with different nucleotide substitution models for each partition.

Bayesian Evolutionary Analysis by Sampling Trees (BEAST 1.8.2; Drummond et al., 2012) was used to estimate the divergence time in Orobanchaceae. Bayesian Evolutionary Analysis Utility (BEAUti 1.8.2) was used to output BEAST input files. Owning to the lack of a reliable fossil record within the family (Soltis et al., 2011; Tank et al., 2015; Uribe-Convers and Tank, 2015), two external fossil calibration points (Call and Dilcher, 1992; Collinson et al., 1993) were used in the dating analysis. Both the stem age of Solanaceae and the age of the most recent common ancestor of Pedicularis and Olea were constrained with the same lognormal distribution prior (an offset of 44.3 Mya, a mean of 1.5, and a standard deviation of 0.5; Zanne et al., 2014). These two nodes are relatively close to Orobanchaceae (Zanne et al., 2014). Ten independent Monte chains Carlo Markov (MCMC) were conducted to ensure convergence in divergence time. Each run consisted of 25 million generations (sampling every 1,000 steps) with the GTR + G model, a Yule tree prior and an uncorrelated lognormal clock. Tracer 1.6.0 (Rambaut and Drummond, 2012) and an R package, RWTY (Warren et al., 2017), were used to assess convergence and stationarity of each MCMC chain. The RWTY test showed that each MCMC chain reached convergence within 5 million generations. Therefore, our analyses of 25 million generations guarantee convergence and stationarity of each MCMC chain. The samples from each run were combined by LogCombiner 1.8.2, until ESS ≥ 200. The final maximum clade credibility tree was generated using TreeAnnotator 1.8.2 after removing 10% of the samples as burn-in. The topology and node height with 95% highest posterior density (HPD) were visualized in FigTree 1.3.1 (Rambaut, 2006).

Of the genes used in our ML analyses, the nuclear gene ITS has a limited capacity to resolve the phylogeny of Orobanchaceae at the genus level (Schneeweiss et al., 2004a; Wolfe et al., 2005; Gussarova et al., 2008). Our ML tree based on ITS data shows good support in four clades (I, IV, V, and VI), but poor support for Clade II [maximum likelihood bootstrap (MLBS) < 50%]. Although Clade III collapses in the ML tree, section Orobanche forms a clade, in which the relationships among most species obtain strong support (Figure S1). The plastid gene rps2, which is known to be retained in all hemi- and holoparasites in Orobanchaceae, behaves similarly to ITS. The coding gene matK has been widely used to infer phylogenetic relationships among closely related groups. Compared with rps2, combined plastid data (rps2 + matK) can resolve the basal clades and support the relationships of the species within the main clades well. The nuclear genes PHYA and PHYB are associated with seed germination (Mathews and Donoghue, 2000; Franklin and Whitelam, 2005). Previous studies suggested that PHYA is the most useful single locus to resolve the phylogenetic relationships within Orobanchaceae (Mathews and Donoghue, 2000; Bennett and Mathews, 2006), but the phylogenetic trees of PHYB, by contrast, showed higher resolution of the relationships among and within the main clades (Figures S5, S6, Table S4). In general, phylogenetic resolution of the combined data is superior to that of any single gene. To gain a better understanding of the phylogenetic relationships among the species of Orobanchaceae, it is necessary to combine nuclear genes with plastid genes as in the previous study by McNeal et al. (2013). In this study, our phylogenetic analyses based on the combined data show strong support among and within the major clades in Orobanchaceae (Figure 2).

All ML trees based on single or combined genes support that P. tubiflorus clusters with the members of the Orobancheae (Figures 1, 2, Figures S1–S11). Traditionally, the genus Orobanche was divided into four sections: Gymnocaulis, Myzorrhiza, Trionychon, and Orobanche. As far as Orobanche in Clade III is concerned, our phylogenetic analyses based on both single and combined genes obtained the same results as the previous study by Schneeweiss et al. (2004a), i.e., Orobanche falls into two lineages: the subgenus Orobanche which contains the sections Orobanche and Diphelypaea, and the subgenus Phelipanche which contains the sections Gymnocaulis, Myzorrhiza, and Trionychon. Almost all ML trees in this study strongly support P. tubiflorus as the closest relative of section Orobanche (99 or 100% MLBS), with the exception of the ITS tree (72% MLBS). Independent analysis of either ITS or rps2 supports that P. tubiflorus is nested in section Orobanche (Figure 1, Figures S7–S10). Although we looked for morphological differences between Phacellanthus and the section Orobanche, we found only one: the dehiscing capsule has three valves in Phacellanthus but two valves in the section Orobanche. Except for this difference, they present high similarities in the shape of the stem, leaf, calyx, corolla, and capsule, the arrangement of the leaves, the number of stamina, and so on. Considering all the evidence above, it is better to merge Phacellanthus into section Orobanche.

Unlike the situation of Phacellanthus, almost all analyses from single and combined genes support G. mupinensis as an independent lineage. Even if ITS shows poor resolution in Orobancheae, G. mupinensis forms a polytomy with several large clades (Figure S1, Table S4). ML trees constructed with PHYA, PHYB, two plastid genes or combined data from five genes (Figures S5, S6, Table S4) group G. mupinensis within Clade III, near to Orobanche, but G. mupinensis is not nested in the clade including Boschniakia, Epifagus, and Conopholis. The phylogenetic positions of Gleadovia and Phacellanthus clearly support the tribe Orobancheae revised by McNeal et al. (2013), but do not support the tribe Gleadovieae which groups Gleadovia with the genera Mannagettaea, Phacellanthus, and Christisonia (Zhang and Tzvelev, 1998). Indeed, only one chamber is fertile in Christisonia and Aeginetia, supporting their sister-group relationship in Clade VI. However, the anthers of the species of Gleadovia and the members of tribe Orobancheae are fertile and have two equal chambers. By contrast, Gleadovia shares more morphological similarities with Mannagettaea and Phacellanthus despite the slight differences in the number of placentae among these three genera and the number of carpels between Gleadovia and Phacellanthus.

Several species in Orobanchaceae such as A. orobanchoides, S. gesnerioides, S. hermonthica, and T. alpina were regarded as holoparasites before McNeal et al. (2013) showed that these species retain functional chlorophyll and photosynthesize at least in part of their life cycle, although sometimes the photosynthesis rates are very low. Without considering these species, our results support three origins of holoparasitism in Orobanchaceae. We note that certain hemiparasitic genera cannot be included in this study due to the limitation of sampling. Because Clade III is strictly restricted at the circumscription of Orobanchaceae sensu stricto containing only holoparasites, and the holoparasites in Clade V consist of only the species of Lathraea, the holoparasitic clade (Hyobanche, (Harveya, (Aeginetia, Christisonia))) within Clade VI becomes the most likely lineage in which hemiparasitic species might be nested. However, these four closely related holoparasitic genera are morphologically so similar that their holoparasitism is most likely homologous. Any hemiparasites nested among them would therefore have regained the ability to photosynthesize. We consider this highly unlikely because no such reversal from holo- to hemiparasitism has been reported to date; but even in that case, the number of origins of holoparasitism would remain unchanged. Thus, if we accept that A. orobanchoides, S. gesnerioides, S. hermonthica, and T. alpina are hemiparasites (McNeal et al., 2013) and that Eremitilla and Platypholis belong to Clade III, there were three independent origins of holoparasitism in Orobanchaceae, each time from hemiparasites.

Orobanchaceae was estimated to have a mid-Tertiary Laurasian origin (Wolfe et al., 2005; Soltis et al., 2011; Tank et al., 2015; Uribe-Convers and Tank, 2015). Based on the fossil record of close relatives of Orobanchaceae, we can compare the relative divergence times of the holoparasitic clades (Figure 3, Figure S15). Our data show that the oldest holoparasitic clade is Orobancheae. The second holoparasitic clade is (Hyobanche, (Harveya, (Aeginetia, Christisonia))) in Buchnereae. The crown group of this clade seems to postdate those of Orobanche and the clade (Boschniakia, (Kopsiopsis, (Epifagus, Conopholis))). The youngest holoparasitic clade comprises only the genus Lathraea. The age of this genus is almost equal to the divergence between Cistanche phelypaea and C. tubulosa.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.