The substantial morphological and genetic change in parasitic plants compared to their closest relatives has traditionally hindered the inference of evolutionary relationships for several reasons, including: limited comparisons within lineages due to a paucity of characters following “parasite reduction syndrome” (Colwell 1994), difficulty in identifying sister groups due to fewer characters, and for characters that do exist, more extreme variation in character states.

One such clade of parasitic flowering plants is the genus Arceuthobium M.Bieb. Commonly referred to as dwarf mistletoes, they are ecologically and economically important obligate parasites that infect young stem tissue of coniferous hosts (specifically Pinaceae and Cupressaceae) throughout North and Central America, Africa, Asia, and Europe (Hawksworth and Wiens 1996). Dwarf mistletoes obtain the majority of their carbohydrates, mineral nutrients, and water from infected hosts, thereby negatively influencing host carbon and water balances (Hull and Leonard 1964a, 1964b; Hawksworth and Wiens 1996). Severe infestations of dwarf mistletoes can significantly alter forest structure and composition via reduced host growth and longevity as well as increased susceptibility of hosts to insect infestation and catastrophic wildfire (Hawksworth and Wiens 1996; Mathiasen et al. 2008).

From an economic perspective, dwarf mistletoe infections incite growth loss and deform commercially valuable timber species, particularly in the western US. (Hawksworth and Wiens 1996; Mathiasen et al. 2008). Although long considered pests, the conservation and ecological importance of dwarf mistletoes, including influences on stand- to ecosystem-scale dynamics, wildlife diversity, and host habit, has gained considerable attention over the last two and a half decades (Hawksworth and Wiens 1996; Mathiasen 1996; Shaw et al. 2004; Mathiasen et al. 2008).

Despite the uncontroversial generic circumscription of Arceuthobium, three factors (which mirror the general difficulties listed above) contribute to taxonomic controversy at the species level: limited variability in gross morphology, species identifications based on combinations of partially overlapping continuous characters, and limited phylogenetic evidence. Recent debate about dwarf mistletoe taxonomy, particularly in western North America, evinces a deeper dispute regarding the evolution and relationships of these plants. The most recent generic monograph by Hawksworth and Wiens (1996) recognized 42 species and four additional subspecies that were “geographically restricted” and “distinguished by a few relatively small but consistent variations.” Nonetheless, Kuijt (2012) lumped into synonymy eleven taxa from Arceuthobium section Campylopoda in his floristic treatment for California (The Jepson Manual, ed. 2), recognizing only three species: A. americanum Engelm., A. campylopodum Engelm., and A. douglasii Engelm. In response, Mathiasen and Kenaley (2016) provided a summary of the morphological, physiological, and/or host specialization differences that in their view support continued taxonomic recognition of Californian Arceuthobium, as well as its history of classification in western North America. Generally, this evidence is described in the monograph by Hawksworth and Wiens (1996) and subsequent narrowly-focused studies that use univariate and/or multivariate analyses of many morphologic traits to assess inter- and intraspecific boundaries among closely-related taxa (Wass and Mathiasen 2003; Mathiasen and Daugherty 2007, 2009; Scott and Mathiasen 2009; Mathiasen and Kenaley 2015, 2017, 2019; Kenaley et al. 2016). In a few cases, morphometric analyses of Arceuthobium taxa have been supplemented by genetic data (Kenaley and Mathiasen 2013; Reif et al. 2015). Analysis of morphological variation among North American Arceuthobium, particularly controversial members of A. sect. Campylopoda, using ordination and discriminant techniques by Mathiasen and Kenaley (2017, 2019) continues, and has led to the description of additional taxa (Kenaley 2020).

Nickrent (2012, 2016) has taken an intermediate taxonomic position between the aforementioned perspectives by recognizing most Arceuthobium sect. Campylopoda taxa at the sub-specific rank. This decision was largely informed by his interpretation of the available biosystematic and phytogeographic data in the context of the most comprehensive molecular phylogenetic study to date. In that study, Nickrent et al. (2004) sampled the internal transcribed spacer (ITS) as well as the plastid region encompassing tRNA genes trnT, trnL, and trnF (henceforth, trnT-L-F) for each of the 42 species recognized by Hawksworth and Wiens (1996). Their result was a moderately resolved phylogeny and updated sectional classification. Though infraspecific sampling was not sufficient to resolve the monophyly of minimally ranked taxa, Nickrent et al. (2004) found strong support for many infrageneric clades and were able to erect a robust natural, sectional classification. At the same time, closely related taxa (e.g. species within A. sect. Campylopoda) were found to have nearly identical sequence data at the ITS and trnT-L-F loci resulting in a large polytomy. Along with morphological, host, and geographical similarity, some investigators cite this as evidence that the number of Arceuthobium species in western North America, or even taxa altogether, has been overstated (Kuijt 2012; Nickrent 2012).

Equally important in resolving evolutionary relationships, the proper rooting of Arceuthobium among the Viscaceae remains unclear. Nickrent and García (2009) resolved sister Old World and New World clades in ITS and cpDNA analyses rooted by Arabidopsis, yet perplexingly, an outgroup was never used by Nickrent et al. (2004) in their earlier, more comprehensive study of Arceuthobium species relationships. These authors ambiguously stated in the Methods that they used “Old World taxa… to root phylogenetic trees” based on an earlier study of relationships in Viscaceae using 18S rDNA that “placed A. oxycedri basal and sister to A. verticilliflorum and A. pendens.” However, Having only sampled one individual of each species (one from the Old World and two from North America), it would be possible to conclude that all Old World taxa are monophyletic or sister to the New World taxa. Compounding the confusion, the designated root does not appear on the figures presented by Nickrent et al. (2004) or associated text, yet clades and sister taxa, as well as the polarity of character state changes, can only be deduced from rooted trees (Wilkinson et al. 2007).

Therefore, the objective of this study was to infer the most comprehensive phylogeny of Arceuthobium to date. Our specific goals were to 1) increase sampling, emphasizing multiple individuals at the infraspecific level and within Arceuthobium sect. Campylopoda from distinct hosts or geographic ranges in order to test monophyly of species, 2) increase the amount of sequence data from both the nuclear and plastid genomes for each sample to more accurately estimate topology, branch lengths, and infraspecific relationships, and 3) correctly root the genus using plastome and nrDNA sequence data.

Materials and Methods

Taxonomic Sampling and DNA Extraction—Seventy specimens from the genus Arceuthobium representing all species recognized by Hawksworth and Wiens (1996) or subsequent authors were selected for DNA sequencing, if possible from distinct hosts, geographic localities, and infraspecific taxa. In some cases, DNA aliquots from previous studies were available (e.g. Reif et al. 2015), but in most cases tissue was sampled from herbarium specimens and ground using mortar and pestle, taking care to avoid flowers or fruits if possible (Appendix 1).

DNA extractions was extracted using a DNeasy Plant Mini Kit (Qiagen, Valencia, California) with slight modifications to manufacturer's instructions. Eluted DNA was quantified using a Nanodrop spectrophotometer or Qubit 3.0 fluorometer (dsDNA assay; Thermo Fisher Scientific, Waltham, Massachusetts). Raw extractions were diluted or concentrated to 1.3 µg DNA in 100 µL solvent. In some cases, as little as 0.7 µg was used due to low extraction yield. The DNA was sheared to 300–500 bp fragments using a Bioruptor sonicator (Diagenode, Denville, New Jersey) at low power for 5–10 minutes of 30s on/30s off cycles. Samples were electrophoresed and visualized on an agarose gel using ethidium bromide to verify the approximate size of DNA fragments. Genomic DNA (gDNA) libraries were prepared in-house using a modified version of Meyer and Kircher's (2010) protocol. Thereafter, samples were each sequenced in 1/50 to 1/40 of an Illumina HiSeq 4000 DNA sequencer lane (2 × 150 bp paired end format) at the Vincent J. Coates Genomics Sequencing Laboratory (Berkeley, California). The gDNA libraries had an average size of 440 bp (Bioanalyzer, Agilent, Santa Clara, California) and were multiplexed along with samples from other projects across two lanes. Due to a small number of samples that failed the library prep or sequencing phases, four species are only represented in this study by preexisting data from Genbank: Arceuthobium azoricum Hawksw. & Wiens, A. divaricatum Engelm., A. pusillum Peck, and A. tibetense H.S.Kiu & W.Ren. Additionally, Kenaley (2020) described two new subspecies of A. abietinum (Engelm.) Hawksw. & Wiens after sequencing had been completed, thereby precluding complete infraspecific sampling.

DNA Sequence Assembly and Alignment—Sequence reads were demultiplexed by the sequencing center. Genomic regions selected for analysis included most of the nuclear ribosomal cistron (ETS, 18S and 5.8S rDNA, both ITS regions, and a small portion of 26S rDNA) and several plastome regions including both coding and noncoding sequences (Table 1) assembled using GetOrganelle v. 1.7.0 (Jin et al. 2018). Though we were generally able to assemble the full nrDNA cistron, the 28S rDNA gene and most of the intergenic spacer were omitted from phylogenetic analysis due to minimal variation and difficulty in confident alignment, respectively.

Resultant sequences were consistent with, but superior in quality to, preliminary reference-guided assembly results. In this pipeline, we quality-filtered the raw reads using Trimmomatic v. 0.36 (Bolger et al. 2014) to identify seed matches with the index adapters (16 bases) allowing maximally one mismatch. These seeds were extended and clipped if in the case of paired end reads a score of 30 was reached (∼50 bases), or in the case of single ended reads a score of 10, (∼17 bases). Only paired reads greater than 40 bp with an average quality per base of at least 20 in a 10 bp sliding window were retained. Assembly to within-genus references was then performed using BWA v. 0.7.17 (Li and Durbin 2009) and subsequent filtering using Samtools v. 1.9 (Li et al. 2009) based on a Phred quality of 20 and read depth of at least 15.

Alignments for each gene were made using MUSCLE v. 3.7 implemented either in Geneious or the CIPRES portal, then visually inspected and refined if necessary. Outgroup taxa were added from preexisting Viscaceae sequences available on GenBank, specifically Viscum album L., Phoradendron leucarpum (Raf.) Reveal & M.C.Johnston, Osyris alba L., and Osyris lanceolata Hochst. & Steud. ex A.DC. (Genbank accession numbers appended to tip labels on Figs. 1–3.). High levels of sequence divergence hindered homology assessment across all sites. Therefore, only the more conserved nrDNA regions (5.8S, 18S and part of the 26S) and open reading frames of the plastome genes petA, rbcL, accD, and ycf2 from outgroup taxa were aligned with the complete sequences of ingroup samples. The total length of our genus-wide concatenated alignment was 38,871 bp, of which 3003 bp were from the nuclear ribosomal cistron and 35,868 bp from several plastome regions (Table 1).

To better infer relationships among the complex of taxonomically controversial species in Arceuthobium sect. Campylopoda, we assembled and aligned a 30,159 bp region of the large single-copy region using GetOrganelle (Jin et al. 2018) and MUSCLE v. 3.7 implemented on the CIPRES portal (Plastome contig 6 in Table 1). A region that aligned poorly, comprising 1140 positions between the genes atpB and rps4, was removed prior to analysis. To the remaining 29,016 bp, we concatenated two additional contigs from the inverted repeat (Plastome contigs 4 and 5 in Table 1) for a total of 44,659 bp of plastome sequence. This is about half the total unique length of the plastome: For comparison, the large and small single copy regions plus a single inverted repeat of Arceuthobium sichuanense (H.S.Kiu) Hawksw. & Wiens is 86,481 bp (Chen et al. 2020).

Collectively, five different configurations of these data were used in phylogenetic analysis (Table 1): 1) Concatenations of five plastome contigs from the large single-copy and inverted repeat regions (cpDNA); 2) a larger plastome alignment for samples from Arceuthobium sect. Campylopoda (Campylopoda cpDNA); 3) most of the nuclear ribosomal cistron (nrDNA); 4) our nrDNA dataset aligned with all ITS sequences generated by Nickrent et al. (2004) (nrDNA+Nickrent); and, 5) a concatenated nuclear ribosomal cistron and chloroplast genome analysis for A. sect. Campylopoda (Campylopoda nr+cpDNA). Concatenation scheme 4 (nrDNA+Nickrent) was included in order to provide a cladogram with maximal taxonomic sampling at the species and population levels. In this particular analysis, branch lengths may be distorted due to differences in sequencing methodologies (Illumina vs. Sanger) and total sequence length between this study and the ITS-only nrDNA data of Nickrent et al. (2004).

Table 1.

Genomic regions used for phylogenetic inference. Alignment files and supplementary analysis output available on Dryad (Schneider et al. 2021). LSC = Large single copy region. ¹Bold denotes sense strand, plain text with' denotes antisense strand, < or > denotes partial gene, and ^ψ denotes a pseudogene. ²1140 positions between the genes atpB and rps4 were difficult to confidently align, therefore this region was omitted from phylogenetic analysis.

Phylogenetic Analysis—The most appropriate model of nucleotide substitution and rate-heterogeneity that could be implemented in RAxML was GTR + G for the two analyses of Arceuthobium sect. Campylopoda, and GTR + I + G for all others. This was determined by a small sample size corrected AIC model selection criterion when running a greedy search using PartitionFinder2 v. 2.1.1 on XSEDE (Lanfear et al. 2016) through the CIPRES portal. PartitionFinder2 also identified the best gene region partitioning scheme for our nrDNA and combined datasets. Maximum likelihood trees for each data set were inferred using RAxML-HPC v. 8.2.12 on XSEDE (Stamatakis 2014) using the “CAT” approximation for analyses with > 50 tips as recommended by the user manual. Rapid bootstrap analysis was performed for 1000 iterations.

Data Access—The following data and analysis files are freely available on Dryad (Schneider et al. 2021): 1) all DNA sequence alignments; 2) partitioning schemes and models of nucleotide substitution and rate-heterogeneity used for each analysis; 3) maximum likelihood phylograms for analyses shown in Figs. 1–3 (Newick format); 4) results of concatenated Campylopoda nr+cpDNA analysis (Newick and .pdf format); and, 5) voucher information for all 70 successfully sequenced samples (Appendix 1) in alternative formats (.csv, .xls, and .pdf). DNA aliquots are archived with the fifth author (SCK) and available upon request.

Results

Congruent Plastome and Nuclear Support for Most Infrageneric Relationships—Both plastome (cpDNA) and nrDNA data strongly supported the monophyly of Arceuthobium and its two subgenera, the Old World Arceuthobium subgenus Arceuthobium and the New World Arceuthobium subgenus Vaginata Hawksw. & Wiens (Figs. 1–3, BS = 100%). Across all analyses, the clades representing each subgenus had relatively long stem branches. These branches accounted for approximately 1/3 to 2/3 of the molecular evolution since the most recent common ancestor of Arceuthobium, depending on the species and genomic compartment (plastid or nuclear).

All eleven of the taxonomic sections proposed by Nickrent et al. (2004) were strongly supported as monophyletic based on nrDNA data (BS ≥ 86, Figs. 2–3). Five of these sections, A. sect. Arceuthobium, A. sect. Americana, A. sect. Campylopoda, A. sect. Globosa, and A. sect. Penda, were also resolved as monophyletic with strong support using plastome data (BS = 100, Fig. 1A). By contrast, plastome relationships strongly supported the non-monophyly of A. sect. Rubra and A. sect. Vaginata (Fig. 1A, see section below). Plastome monophyly could not be assessed for the remaining three sections (A. sect. Azorica, A. sect. Minuta, or A. sect. Pusilla) due to limited sampling.

Several relationships among the taxonomic sections were supported by both plastome and nrDNA partitions. Arceuthobium section Americana was resolved sister to a clade comprising all other New World sections (BS = 100%, Figs. 1–3). Likewise, A. sect. Rubra and A. sect. Vaginata were resolved as monophyletic from nrDNA data (BS ≥ 85%, Figs. 2–3), and in a polytomy with A. sect. Pusilla when using plastome data (Fig. 1A). However, relationships among several New World sections could not be resolved by nrDNA data. Arceuthobium section Campylopoda, A. sect. Globosa, A. sect. Pusilla, and the clade of A. sect. Rubra + A. sect. Vaginata formed a large polytomy (BS ≤ 60%, Figs. 2–3), though plastome data supported a clade comprising A. sect. Pusilla, A. sect. Rubra, and A. sect. Vaginata (BS = 100%).

We were unable to generate any usable sequence data from only one section, the monotypic Arceuthobium sect. Azorica. However, adding existing ITS sequences from this section and others to our nrDNA analysis provided modest support that all other Old World taxa (i.e. A. sect. Chinense and A. sect. Arceuthobium) form a clade (BS = 75%) sister to the narrow endemic A. azoricum (Fig. 4).

Congruent Species Relationships—Many species relationships within sections were resolved with strong support (BS = 100% unless otherwise noted; Figs. 1–3, summarized in Fig. 4). Among Old World taxa and within A. sect. Arceuthobium, samples of A. juniperi-proceri Choiv. were sister to two samples of A. oxycedri (DC.) M.Bieb. based on our nrDNA data (Fig. 2). However, plastome data did not provide appreciable support for these relationships (Fig. 1A). Despite being represented only by the ITS sequence of a single specimen, Arceuthobium tibetense was resolved in a clade with, but indistinct from, most other samples of A. sect. Arceuthobium (Fig. 3, BS = 90%). We found strong support for a monophyletic A. pini Hawksw. & Wiens with the expanded nrDNA+ Nickrent dataset (Fig. 3, BS = 88%) and the more inclusive clades adding monophyletic species A. sichuanense and A. chinense Lecomte, respectively (Figs. 1–3). Samples of A. minutissimum Hook.f. also were monophyletic (Fig. 3, BS = 100%), although the placement of this long-branch species within A. subgenus Arceuthobium was uncertain (Fig. 4; see below).

Fig. 1.

Maximum likelihood phylograms of Arceuthobium inferred from plastome data. Bootstrap values ≥ 70 are shown. A. Genus-wide analysis showing relationship among species and taxonomic sections (sidebars at right). Branch lengths connecting outgroup taxa (in gray) are drawn arbitrarily, all others to scale. B. Separate analysis of A. sect. Campylopoda. Host affinities shown for specimens of taxa not resolved as monophyletic, all others in Appendix 1. Outgroup (A. blumeri) not shown. Scale bar represents 0.001 substitutions per site, or approximately 29 substitutions across the alignment.

Among New World species, analysis of both plastome and nrDNA data supported A. abietis-religiosae Heil as sister to the other two species of A. sect. Americana, A. verticilliflorum Engelm. and A. americanum (nrDNA BS ≥ 84%; plastome BS = 100%), and all three of these species were supported as monophyletic (BS = 100%, Figs. 1–3). Within A. sect. Rubra, A. yecorense Hawksw. & Wiens resolved sister to A. oaxacanum Hawksw. & Wiens + A. rubrum Hawksw. & Wiens. Arceuthobium gillii Hawksw. & Wiens and A. nigrum (Hawksw. & Wiens) Hawksw. & Wiens resolved sister in all analyses (BS ≥ 95%, Figs. 1–3). In the closely-related A. sect. Vaginata, A. hondurense Hawksw. & Wiens was sister to a clade of A. vaginatum (Humb. & Bonpl. ex Willd.) J.Presl samples (BS = 100% Figs. 1–3).

Fig. 2.

Maximum likelihood phylogram of Arceuthobium inferred from the nuclear ribosomal cistron (ETS–18S–ITS). Bootstrap values > 70% shown. Branch lengths connecting outgroup taxa (in gray) are arbitrary; all others to scale. Taxonomic sections indicated by sidebars at right.

Fig. 3.

Maximum likelihood phylogram of Arceuthobium inferred from the nuclear ribosomal cistron. Bootstrap support values ≥ 75% shown. Tip labels in black font represent sequences newly generated for this study (see Fig. 2). Tip labels in gray font represent ITS sequences generated by Nickrent et al. (2004), with Genbank accession number indicated after the taxon name. Asterisk indicates subspecific determinations made by the last author of this study (RLM). Branch lengths connecting outgroup taxa (also in gray) are arbitrary; all others to scale.

Fig. 4.

Summary cladogram of Arceuthobium based on plastome and nrDNA sequence data. Nodes with BS < 80% (nrDNA) or BS < 90% (cpDNA) collapsed into polytomies, except as marked. Plastome and nrDNA data support conflicting relationships regarding the position of A. minutissimum and some species within the clade of New World taxa shaded in gray, but elsewhere are congruent. Asterisk indicates four taxa represented only by ITS sequences published by Nickrent et al. (2004). For relationships within A. sect. Campylopoda, see Fig. 1B.

Seven other species were resolved as monophyletic in both nrDNA and plastome analyses, namely A. globosum Hawksw. & Wiens, A. strictum Hawksw. & Wiens, A. gillii, A. hondurense, A. vaginatum, A. blumeri A.Nelson, and A. microcarpum (Engelm.) Hawksw. & Wiens (BS = 89% for A. microcarpum nrDNA, otherwise 100%; Figs. 1–2). Samples from each of the species A. bicarinatum Urban, A. divaricatum, A. douglasii, A. oaxacanum, and A. pusillum. were resolved as monophyletic in the expanded nrDNA+Nickrent analysis (BS ≥ 98%, Fig. 3), but could not be confirmed with plastome data due to insufficient sampling. Conversely, plastome data provided strong support for the monophyly of A. nigrum and three species from A. sect. Campylopoda (A. apachecum Hawksw. & Wiens, A. cyanocarpum (A.Nelson ex Rydb.) Coulter & A.Nelson, and A. microcarpum; BS ≥ 99%, Fig. 1B), but insufficiently resolved by nrDNA data (Figs. 2–3). Plastome sequences of all A. tsugense (Rosen.) G.N.Jones samples except the one identified as subsp. mertensianae Hawksw. & Nickrent formed a strongly supported clade along with one sample of A. laricis (Piper) St.John (BS = 100%, Fig. 1B), but not in the nrDNA analyses (Figs. 2–3). The sample of A. laricis, from Shoshone County, Idaho, was collected on a secondary host, Picea engelmannii, “near heavily infected mountain hemlock and western larch [Tsuga mertensiana and Larix occidentalis]” (Mathiasen 9302; UC 2056905).

Infraspecific Relationships—The nrDNA+Nickrent analysis also resolved A. rubrum as paraphyletic with respect to A. oaxacanum (BS = 100%, Fig. 3); however, plastome sampling was insufficient to confirm this result. Within A. sect. Campylopoda, all taxonomic species not mentioned above were not supported as monophyletic or could not be assessed due to only a single sample (Fig. 1B). The plastome analysis did resolve a clade comprising all samples of A. monticola Hawksw., Wiens & Nickrent plus one sample of A. siskiyouense Hawksw., Wiens & Nickrent (RLM 05–21, BS = 100%, Fig. 1B), but this clade was not supported by nrDNA data (Figs. 2–3).

The monophyly and relationships of infraspecific taxa were less clear. If monophyletic (nrDNA BS = 87%, Fig. 3), Arceuthobium hondurense Hawksw. & Wiens subsp. hawksworthii (Wiens & C.G.Shaw bis) Mathiasen is nested within samples of A. hondurense subsp. hondurense (BS ≥ 98%, Figs. 1A, 3) but with minimal sequence differentiation (i.e. branch lengths). Our four samples of A. globosum showed strongly supported phylogenetic structure in nrDNA and plastome analyses, but the inferred subspecific relationships conflicted between the two datasets (BS ≥ 81%). Moreover, when the ITS-only samples of Nickrent et al. (2004) were included, only A. globosum subsp. globosum was strongly supported as monophyletic (BS = 96%, Fig. 3). Within A. vaginatum, the clade representing A. vaginatum subsp. durangense Hawksw. & Wiens was nested within a well-supported grade of two A. vaginatum samples in the plastome analysis (Fig. 1); however, nrDNA patterns were more complex (Figs. 2–3). While samples of both species collectively formed a strongly supported clade (BS = 100%), any modestly supported substructure was not clearly associated with taxonomic assignments. Three samples of A. vaginaum subsp. durangense formed one subclade (BS = 81%), while two others, including the isotype, resolved within a separate clade with an ITS clone from A. vaginatum subsp. vaginatum (Genbank AY288286) (BS = 93%, Fig. 3).

Within Arceuthobium sect. Campylopoda, plastome data supported a monophyletic subspecies A. microcarpum subsp. microcarpum sister to a single sample of A. microcarpum subsp. aristatae J.M.Scott & Mathiasen (BS = 100%, Fig. 1B), but our three samples were identical at the nrDNA locus, with 100% pairwise identity over all 2836 sites. Also notable, our three samples of Arceuthobium abietinum did not form a clade. One sample each of A. abietinum subsp. wiensii Mathiasen & C.M.Daugherty and A. abietinum subsp. magnificae Mathiasen & Kenaley formed a strongly supported clade (BS = 100%, Fig. 1B) but the sample of the newly described A. abietinum subsp. mathiasenii Kenaley was resolved sister to A. microcarpum (BS = 99%). Arceuthobium abietinum subsp. mathiasenii primarily parasitizes Abies, but in Mexico occasionally infects Picea mexicana Martinez and rarely Pinus strobiformis Engelm. Like nearly all taxa within A. sect. Campylopoda, nrDNA sequences of A. abietinum were indistinct (Fig. 2).

Three Instances of Plastome and nrDNA Incongruence—Two instances of discordance between relationships inferred from plastome and nrDNA were well-supported, and discordance in the position of A. minutissimum was modestly supported (Fig. 4). Both cases of well-supported discordance were within the New World clade (Arceuthobium subgenus Vaginata). First, nrDNA data supported A. sect. Penda as sister to all New World sections except A. sect. Americana (BS = 100%, Figs. 2–3), whereas plastome data supported a clade of A. sect. Penda, A. sect. Rubra, A. sect. Vaginata, and A. sect. Globosa sister to a clade comprising A. sect. Minuta and A. sect. Campylopoda (BS = 100%, Fig. 1). Second, the monophyly of A. sect. Vaginata and A. sect. Rubra and their sister relationship was supported by nrDNA (BS ≥ 85%, Figs. 2–3), but plastome data was more complicated. From A. sect. Vaginata, A. vaginatum was resolved sister to three species from A. sect. Rubra (A. rubrum, A. oaxacanum, and A. yecorense), whereas the two remaining species in A. sect. Rubra (A. gillii and A. nigrum) were nested within a grade of the remaining species in A. sect. Vaginata (A. hondurense and A. strictum, respectively; BS ≥ 85%, Fig. 1).

Finally, the position of the Old World species A. minutissimum remained ambiguous. The nrDNA analysis resolved this species sister to the clade A chinense, + A. sichuanense + A. pini, thereby supporting a monophyletic A. sect. Chinense (BS ≥ 86%, Figs. 2–4). Plastome data provided modest support for A. minutissimum as sister to A. sect. Arceuthobium (BS = 80%, Figs. 1, 4). Separate by-contig analyses indicated this signal appeared to be limited to the Contig 2 portion of the large single-copy region (Table 1). Contig 4 provided very weak support for the nrDNA hypothesis (BS = 68%) and the three longest contigs (1, 3, and 5) provided no support for any particular relationship among these three groups (BS ≤ 48%, data not shown).

Minimal Pseudogenization of Plastome Genes—While not a main focus of this study, the large plastome contigs revealed uniformity in plastome gene content across Arceuthobium and high similarity to autotrophic relatives. In all samples, only cemA and ndhB were pseudogenized (Table 1). Otherwise, open reading frames for coding genes were intact. Substantial insertions and deletions were uncommon, except for within the 16S-23S rDNA spacer. This region also showed the greatest variability of plastome reduction, ranging from 2033bp in A. gillii to 486–487 bp in A. hondurense, and 982–1049 bp among Old World taxa.

Discussion

Evolution and Biogeography of Arceuthobium—Old World (A. subgenus Arceuthobium) and New World (A. subgenus Vaginata)taxa each form two strongly supported and reciprocally monophyletic lineages based both on plastome (Fig. 1) and nuclear ribosomal data (Figs. 2–3), with the genus rooted at the most recent common ancestor of these two lineages. However, the biogeographical history of stem-lineage Arceuthobium remains unclear due to uncertainty in the identity of its closest extant relatives and conflicting circumstantial evidence discussed below.

The fossil record indicates that diversification of Arceuthobium in Europe had begun by the late Eocene (Lutetian to Priabonian, 48–34 Ma; Sadowski et al. 2017). Many of the fossil Arceuthobium specimens contain characters that were not present in extant species and likely are plesiomorphic, such as expanded, oblanceolate leaves, a unique distribution of stomata, and unfused squamate bracts (Sadowski et al. 2017). These features support a hypothesis that stem-lineage Arceuthobium likely was present in the Old World, and that by the late Miocene plants directly attributable to A. oxycedri had evolved (Łańcucka-Środoniowa 1980). Concordant with an Old World origin of Arceuthobium, the earliest known New World Arceuthobium fossils date only to the Miocene. However, the three known localities of Arceuthobium pollen are scattered across western North America, in the Alaska Range (Wahrhaftig et al. 1969), Colorado (Weber 1965), and Wyoming (Leopold and Denton 1987), indicating that dwarf mistletoes were geographically widespread during the Miocene.

Nickrent and García (2009) advocated a New World origin for Arceuthobium due to extant species diversity, increased plas tome reduction (i.e. deletions in the 16S-23 Sintergenic spacer of Old World species), and greater sequence similarity between New World taxa to outgroups (cf. Results, Minimal Pseudogenization of Plastome Genes, above). However, we find these arguments unconvincing. First, the lack of formal analysis of diversification rates or potential rate shifts aside, the high extant diversity in the Mexican highlands cited by Nickrent and García (2009) could equally be a result of increased rates of recent speciation (i.e. neoendemism) or reduced extinction rates in climate refugia. Second, additional sampling at the 16S-23S spacer revealed both the most intact (i.e. longest) and reduced sequences are in the New World clade; however, across larger regions of the plastome sections levels of plastome reduction were minimal. Moreover, given that New World and Old World Arceuthobium comprise sister lineages (Figs. 1–3), differences in plastome loss must have happened after the lineages diverged and presumably dispersed to their respective hemispheres. This is to say, genomic differences in reciprocally monophyletic clades cannot be evidence for or against a biogeographical hypothesis.

Relationships among extant genera of Viscaceae, many of which are limited to particular hemispheres, remain highly uncertain. While a recent analysis using ITS showed strong support for Arceuthobium sister to the predominantly southeast Asian and Australian genus Korthalsella Tiegh. (BS ≥ 96%, posterior probability (PP) = 1, Fig. S3 of Maul et al. 2019), the most recent multi-locus analyses including some plastome regions indicate that the Old World genera Korthalsella and Ginalloa Korth. represent a clade sister to the predominantly New World Dendrophthora Eichler and Phoradendron Nutt. (PP > 0.9, Nickrent et al. 2019; see “Viscaceae_Amphorogynaceae_ cp.tre” in Schneider et al. 2021). Relationships among other Viscaceae genera Viscum L., Notothixos Oliv., and Arceuthobium also were ambiguous (PP < 0.7). These results conflict with other few-gene analyses using different loci (e.g. combined analyses of Maul et al. 2019, and Der and Nickrent 2008) in which Arceuthobium weakly resolves with the Viscum + Notothixos clade. In short, genomic-scale data is clearly needed to clarify generic relationships in the Viscaceae. Until such data can be obtained, thus allowing hypotheses involving the origin and dispersal of Arceuthobium to be rigorously tested and statistically evaluated, the biogeographical history of dwarf mistletoes remains an open question. Therefore, we find the fossil evidence, however tentative, to be most persuasive at present in support of the hypothesis that Arceuthobium evolved in the Old World.

Sectional Relationships and Plastome/nrDNA Incongruence—Consistent with previous morphological and molecular evidence underlying the sectional taxonomy of Nickrent et al. (2004), we found strong support for monophyletic sections, except among Arceuthobium sect. Rubra and A. sect. Vaginata in the plastome analysis only (Figs. 1–3). Surprisingly, additional data did not aid in resolving several polytomies among New World sections that were first reported by Nickrent and García (2009). However, our results do confirm several instances of nrDNA/plastome incongruence within the New World clade sister to A. sect. Americana (Fig. 4). By contrast, we propose that incongruence in the placement of A. minutissimum (Figs. 1–4) is due to incomplete lineage sorting or insufficient phylogenetic signal. Most regions of the plastome were phylogenetically uninformative at this node or showed conflicting but weakly supported alternative topologies (data not shown).

The alternative positions of Arceuthobium minutissimum represent alternative configurations for rooting the Old World A. subgenus Arceuthobium. Such uncertainty is not surprising on account of long stem branches for both A. minutissimum and A. subgen. Arceuthobium. Our nrDNA analyses (Figs. 2–3), which place A. minutissimum within A. sect. Chinense, are consistent with morphological, ecological, and biogeographical evidence used by previous taxonomists (Hawksworth and Wiens 1996). Arceuthobium minutissimum is found in the Himalayas, geographically intermediate between A. sect. Arceuthobium, found from Pakistan westward through central Asia, the Mediterranean Basin, and eastern Africa (Hawksworth and Wiens 1996), while the rest of A. sect. Chinense is found in the Himalayas and eastward throughout China. Like other taxa in A. sect. Chinense, A. minutissimum parasitizes Pinaceae, whereas members of A. sect. Arceuthobium and A. sect. Azorica parasitize hosts in Cupressaceae.

Species Concepts and Fine-Scale Taxonomic Implications—In the most recent monographs of Arceuthobium, Hawksworth and Wiens (1972, 1996) relied primarily on a morphological species concept, recognizing “population systems that exhibit suites of characteristics that remain constant within prescribed limits of variation, from generation to generation, on different hosts, and when they occur with other taxa.” In the last decades, increasingly available molecular phylogenetic data coupled with theoretical and computational advancements for inferring phylogenetic trees has facilitated the application of a phylogenetic species concept by improving diagnosis of monophyletic groups (De Queiroz 2007). For example, Nickrent et al. (2004: Table 2) proposed a phylogenetic classification of species, based on ITS and trnT-L-F gene trees. The additional sequence data and population sampling we provide in this study supports the monophyly of each of their proposed species. However, our findings have also revealed additional phylogenetic structure in the genus, both among superspecific groups and towards the tips of our phylogeny. Many of these fine-scale entities corresponded to previously recognized taxa and conform to the criteria of “separately evolving metapopulation lineages” as evidenced by various lines of morphological, ecological (i.e. host affinities), and, as provided by this study, phylogenetic evidence (De Queiroz 2007). Combining our more finely sampled and resolved phylogenetic findings with the numerous recently published morphological studies, in the following sections we clarify the status of taxa not recognized at the species level in Nickrent et al. (2004: Table 2). While philosophically distinct, our ultimate result is concordant with Hawksworth and Wiens's (1996, p.146) approach that recognized “geographically restricted populations delimited by a few relatively small but consistent variations” as subspecies. That is, in many, but not all cases we recognize these synonymized taxa at subspecific rank.

A. globosum—Although we found strongly supported internal structure among individual samples of each of the four subspecies of A. globosum (A. globosum subsp. aureum (Hawksw. & Wiens) Mathiasen, A. globosum subsp. globosum, A. globosum subsp. grandicaule Hawksw. & Wiens, and A. globosum subsp. petersonii (Hawksw. & Wiens) Mathiasen), plastome and nrDNA topologies supported conflicting infraspecific relationships (Figs. 1–2). At the ITS locus, for which we were able to sample multiple populations of each subspecies, only A. globosum subsp. globosum was strongly supported as monophyletic (BS = 96, Fig. 3). Nevertheless, differences among these taxa are clear. Three examples include: the flowering time of A. globosum subsp. petersonii versus the other subspecies, discontinuous ranges among A. globosum subsp. grandicaule, A. globosum subsp. globosum, and the others throughout Mexico and Guatemala, and various morphological differences (Mathiasen et al. 2008). Thus, we follow Mathiasen et al. (2008) in recognizing taxa previously treated under A. aureum as subspecies of A. globosum.

A. oaxacanum and A. rubrum—Arceuthobium oaxacanum was initially described after highly disjunct populations (by over 1000 km) similar to A. rubrum were discovered in southern Mexico but with reported differences in plant color, plant size, inflorescence branching angle, length of pistillate spikes, and some other features (Hawksworth and Wiens 1989). However, more recently, this taxon has not been accepted due to morphological and phenological similarity and genetic similarity at the ITS and trnT-L-F loci (Nickrent et al. 2004; Mathiasen et al. 2009). Our findings suggest that these groups may be more genetically distinct than previously realized. Samples of A. oaxacanum and A. rubrum demonstrated only 95% pairwise sequence identity across large portions of the plastome and 98% pairwise identity at the ETS locus, comparable to other recognized taxa in the genus. Including the ITS sequences published by Nickrent et al. (2004) in a phylogenetic analysis indicated a strongly supported clade of near-identical A. oaxacanum samples nested within a grade of two populations of A. rubrum (BS = 100%, Fig. 4). Therefore, these Oaxacan populations likely represent the early stages of allopatric speciation following reproductive isolation from a Durangan A. rubrum-like progenitor. If diagnostic morphological characters can be found, the molecular differentiation and unique (albeit not reciprocally occurring) hosts would likely warrant recognition of Oaxacan populations as a subspecies within A. rubrum. However, at present, we decline to erect a new combination and treat A. oaxacanum as a synonym of A. rubrum.

A. gillii and A. nigrum–We confirmed a reciprocally monophyletic relationship between A. gillii and A. nigrum plastomes (BS =100%, Fig. 1) first suggested by the trnT-L-F analysis of Nickrent et al. (2004). This is consistent with, but better resolved than, our nrDNA data, in which a monophyletic A. gillii forms a soft polytomy with A. nigrum samples (BS = 100%, Fig. 2). Morphology, ecology, and phytogeography support the distinctiveness of these populations, which were first described as subspecies but elevated to species rank following additional morphological analysis (Hawksworth and Wiens 1989; Kenaley and Mathiasen 2013). Although no synapomorphies for A. nigrum appear present at the nrDNA locus, we predict that they may exist at other, unsampled nuclear loci.

Nickrent et al. (2004) concluded that A. gillii was nested within a paraphyletic A. nigrum on account of the position of a sample (Nickrent 2041, SIU), that resolved sister to all other samples of both species, in conflict with its position in the trnT-L-F analysis. Subsequent studies of A. nigrum using the ITS locus have not found any other populations with a similar sequence or that otherwise support the “paraphyletic A. nigrum” hypothesis (Kenaley and Mathiasen 2013; also see Fig. 3). We were not able to evaluate the specimen in question, but the population grows in an area that includes A. nigrum, A. vaginatum, and A. globosum subsp. grandicaule in close proximity, and at the specific locality some plants in the area had characteristics of A. vaginatum (R. Mathiasen pers. obs.). This anomalous specimen aside, all available morphological and ecological studies (Hawksworth and Wiens 1965, 1989; Kenaley and Mathiasen 2013) show distinctiveness consistent with our plastome and nrDNA phylogenetic support of monophyly for both A. gillii and A. nigrum (Figs. 1–2). Therefore, both of these taxa must be taxonomically recognized, most appropriately at species rank.

A. hondurense (Including A. hondurense subsp. hawksworthii)—Our findings confirm those first reported by Nickrent et al. (2004) that minimal molecular differentiation exists among these populations, with A. hondurense subsp. hawksworthii nested within samples of A. hondurense subsp. hondurense (Figs. 1–3). Our two samples of A. hondurense subsp. hawksworthii formed a modestly supported clade (BS = 82%) at the ITS locus and can be distinguished based on subtle morphological differences and a distinct, non-overlapping reproductive phenology (Mathiasen 2007b). Therefore, the evolutionary history of A. hondurense subsp. hawksworthii is likely analogous to that of A. oaxacanum: a recently derived segregate. Arceuthobium hondurense subsp. hawksworthii should be accepted at subspecific rank if recognized at all.

A. vaginatum (Including A. vaginatum subsp. durangense and A. vaginatum subsp. cryptopodum)—Plastome and nrDNA each demonstrated internal structure within the strongly supported A. vaginatum clade but these relationships were not consistent. For example, plastome but not nrDNA data resolved a single clade of A. vaginatum subsp. durangense populations (Figs. 1–2). Our expanded sampling analysis provided further evidence of continued gene flow, in which an ITS clone from A. vaginatum subsp. vaginatum (Nickrent 2059, ILL; Genbank AY288286) was resolved within one of the two clades of A. vaginatum subsp. durangense samples (BS = 93%, Fig. 3). The host ranges of all three A. vaginatum subspecies partially overlap (Mathiasen 2007a) and Hawksworth and Wiens (1965) observed numerous intermediate individuals between A. vaginatum subsp. cryptopodum (Engelm.) Hawksw. & Wiens and A. vaginatum subsp. vaginatum in Chihuahua, Mexico. Given previous evidence of morphological and to some degree geographical separation (Hawksworth and Wiens 1972, 1996; Mathiasen 2007a), for now they are appropriately recognized at subspecific rank.

Arceuthobium Section Campylopoda—Substantially increasing the sequence data from previous analyses has revealed previously unrecognized phylogenetic structure in this clade. Our results support previous studies that recovered a monophyletic A. blumeri sister to the rest of the clade (Figs. 1–3; Nickrent et al. 2004), but provide plastome evidence of monophyly in several other taxa, including A. cyanocarpum, A. apachecum, A. microcarpum subsp. microcarpum, and A. microcarpum (BS ≥ 99%, Fig. 1B). Among these only A. microcarpum is supported by nrDNA data (Fig. 2). Reif et al. (2015) found more AFLP variation within populations (70%) and modestly more among populations (17%) than among species (13%) when comparing three of these species that parasitize white pines, A. apachecum, A. blumeri, and A. cyanocarpum. Populations of each of these three species formed strongly supported clades in our plastome analysis (BS = 100%), but as the most divergent, A. blumeri was likely more responsible for the low among-species variation than the other two taxa.

Populations of several other taxa within Arceuthobium sect. Campylopoda did not form strictly monophyletic groups but nonetheless showed some structure associated with taxonomic assignment, as well as clustering with geographically nearby samples. For example, both samples of A. monticola resolved in a strongly supported clade with one sample of A. siskiyouense (BS = 100%, Fig. 1B). These taxa often grow sympatrically within their respective ranges along the California-Oregon border but on separate hosts. The other sample of A. siskiyouense was in a clade with populations of A. campylopodum, A. occidentale, and one sample of A. laricis, all from California or the Pacific Northwest (BS = 94%, Fig. 1B). Furthermore, in our plastome phylogeny, all samples of A. tsugense except subsp. mertensianae formed a strongly supported clade along with the sample of A. laricis collected from secondary host Picea engelmannii Parry ex Engelm. (BS = 100%, Fig. 1B; see Results for details on this specimen). Some of these samples also formed a modestly supported clade in the nrDNA analysis (BS = 79%, Fig. 2).

A similar example of geographic structure was found in A. abietinum. The two samples identified morphologically as A. abietinum subsp. magnificae and A. abietinum subsp. wiensii resolved as a clade were both collected in California, whereas A. abietinum subsp. mathiasenii was collected in Utah and forms a clade with the three A. microcarpum samples, collected across Arizona (BS = 99%, Fig. 1B). These taxa have not been formally compared morphologically, as the assignment of A. abietinum subsp. mathiasenii to A. abietinum was based primarily on host affinity to firs (except for rare infection of Pinus strobiformis by both taxa, and occasional infection of Picea A.Dietr. by A. abietinum subsp. mathiasenii). Arceuthobium abietinum subsp. mathiasenii and A. microcarpum also differ in flower size and the proportion of three- versus four-petaled flowers (Mathiasen et al. 2017; Kenaley 2020), so we believe that the plastome affinity is most likely due to introgression.

Nonetheless, geography does not exclusively explain plastome relationships. For example, the clade of two A. cyanocarpum samples were collected nearly at opposite ends of the range of A. sect. Campylopoda, El Paso County, Colorado (Mathiasen 09–03, ARIZ404136), and Inyo County, California (Mathiasen 08–17, ASC00094829). In these cases, we would expect that additional, more rapidly evolving nuclear loci would support the distinctiveness and monophyly of certain taxa within A. sect. Campylopoda.

Though our plastome results (Fig. 1B, and nrDNA in the case of A. microcarpum, Fig. 2) are the first to show phylogenetic structure within this complex, taxa within A. sect. Campylopoda have a long history of recognition based diagnosable morphological and host/ecological differences (Hawksworth and Wiens 1972, 1996; Mathiasen 1982; Mathiasen and Kenaley 2015, 2016, 2017, 2019; Kenaley et al. 2016; Kenaley 2020). As pointed out by Nickrent (2012), these morphometric studies generally have not encompassed all members of a monophyletic group within A. sect. Campylopoda, and so may overestimate the distinctiveness of a taxon compared to its closest relative. However, keys distinguishing all but the most recently described taxa are available (Hawksworth and Wiens 1996; Reif et al. 2015; Mathiasen and Kenaley 2016, 2017).

Phylogenetic structure revealed by plastome and nrDNA analysis (Figs. 1B, 2) extends our understanding that certain groups have progressed at different rates along their separate evolutionary trajectories in acquiring properties such as monophyly (e.g. A. cyanocarpum, A. microcarpum, and A. apachecum), genetic differentiation despite possible introgression (e.g. A. monticola and A. tsugense), or morphological diagnosability (De Queiroz 2007). From an evolutionary perspective, our phylogenetic data support Nickrent et al.'s (2004) model of “incipient speciation whereby the earliest stages of reproductive isolation and genetic differentiation are taking place” among host races or metapopulations of A. sect. Campylopoda. However, given that at least some of these taxa form monophyletic groups, recognizing all of the taxa in the clade sister to A. blumeri as a single entity (e.g. Kuijt 2012) is no longer tenable.

Treating many of these entities as subspecies as advocated by Nickrent (2012) would require additional nomenclatural combinations for other entities that are distinctive morphologically and in some cases phylogenetically (e.g. subspecies of A. tsugense, A. abietinum, and A. microcarpum), but would be consistent with the degree of genetic divergence shown by the subspecies of other taxa. Two downsides to this approach are that it masks some relationships among subspecies and it may have practical conservation implications (Haig et al. 2006). Alternatively, the approach of Mathiasen and Kenaley (2016) of treating most of these taxa at species rank is philosophically defensible under the “unified species concept” of de Queiroz (2007), which emphasizes the unique evolutionary trajectory of metapopulations. However, under this approach all recognized taxa should be treated as full species.

In summary, the phylogenetic structure within A. sect. Campylopoda (not limited to, but including monophyly; Figs. 1B, 2) provides additional evidence that these entities are separately evolving. We hope this will stimulate additional multispecies coalescent, multi-locus phylogenetic, population genetic analysis focused on A. sect. Campylopoda, and recommend that any future molecular or morphological studies encompass well-supported clades to avoid comparisons of para- or polyphyletic groups.

Minimal Plastome Reduction in Arceuthobium—We found only minor loss or pseudogenization of plastome regions in Arceuthobium, consistent with patterns seen across Viscaceae and other hemiparasitic Santalales (Chen et al. 2020). Previous comparative studies found accelerated substitution rates in Arceuthobium relative to Phoradendron (Nickrent and García 2009) and relative to autotrophic relatives (Chen et al. 2020). However, unlike other endoparasites that form systemic infections and have the most reduced plastomes of any plant (e.g. Apodanthacae and Rafflesiaceae; Bellot and Renner 2016; Molina et al. 2014), Arceuthobium is only partially endoparasitic and retains some capacity for photosynthesis (Hull and Leonard 1964b; further evidence reviewed by Hawksworth and Wiens 1996).

Within the inverted repeat, the only putatively nonfunctional genes we found in Arceuthobium were trnA and trnI from the 16S-23S intergenic region (lost in all samples), first noted by Nickrent and García (2009). Plastome reduction in the inverted repeat was variable (see Results). Across 44kb of the large single-copy plastome region, only two coding genes were absent or pseudogenized: cemA and ndhB (Table 1). The NADH dehydrogenase (ndh) genes have been lost in Viscaceae and are often absent in fully autotrophic plants (Petersen et al. 2015; Chen et al. 2020). By contrast, cemA is more central to photosynthesis, encoding a protein involved in carbon uptake. Knockout experiments in Chlamydomonas Ehrenb. show that this gene is not essential for photosynthesis, but nonetheless crucial for photosynthetic performance in high-light environments (Rolland et al. 1997). Absence of a functional cemA in a fully photosynthetic species has otherwise only been noted in the semi-aquatic herb Saniculiphyllum guangxiense C.Y.Wu & T.C.Ku (Saxifragaceae, Folk et al. 2020). In other heterotrophic plants, the loss of cemA is concomitant with the loss of other essential photosynthesis genes during the transition to holoparasitism or full mycoheterotrophy (Graham et al. 2017; Wicke and Naumann 2018). In Viscaceae, this gene is reported present but evolving significantly more quickly compared to the autotrophic outgroup Vitis rotundifolia Michx (Chen et al. 2020). However, visual inspection of the A. sichuanense sequence from that study (Genbank MN414160) reveals the same pattern of numerous internal stop codons as we report in this study, and thus we interpret the gene as nonfunctional.