Hybridization in the Fringed Orchids: An Analysis of Species Boundaries in the Face of Gene Flow

: Natural hybridization between closely related species in sympatry is an evolutionary process that is common in orchids. Once seen as a threat to parent species, interspeciﬁc genetic change is increasingly viewed as a source of novel variation in some ecological contexts. Terrestrial fringed orchids in the genus Platanthera contain several clades with high genetic compatibility among species and many putative hybrids. We used biallelic SNPs generated with 3RAD sequencing to study the hybrid complex formed from the parent species P. blephariglottis, P. ciliaris, and P. cristata with high resolution. The genetic structure and phylogenetic relationship of the hybrid complex revealed site-dependent gene ﬂow between species. We documented extensive hybridization and cryptic hybrids in sympatric sites. Interspeciﬁc genetic exchange is particularly common between P. blephariglottis and P. ciliaris , with cryptic hybrids among putative P. ciliaris samples being more common than parental assignments in sympatric sites. Hybridization across the triad species complex can reticulate lineages and introduce adaptive alleles. Conversely, it can reduce diversiﬁcation rates and introduce maladaptive alleles. Investigation into whether anthropogenic forces are eroding species boundaries, particularly the permeable P. blephariglottis and P. ciliaris boundary, is appropriate for conservation efforts.


Introduction
Orchids constitute one of the most species-rich plant families, with 25,000-30,000 documented species [1,2].The exact number of orchid species and definitions of species boundaries are sources of contention between so-called lumpers and splitters engaged in taxonomical debates [3].Classical definitions for a biological species center around reproductive isolation [4].Further deliberation recognized that species could remain distinct despite some gene leakage due to introgression, defined as the transfer of alleles from one taxon to another divergent taxon [3,5].This concept of semipermeable species lines applies to several orchid species that appear to have diverged despite ongoing gene flow [3,6].When defining lineages as separate species, the extent to which genetic exchange can occur, especially between recently diverged populations, is controversial.
Orchids have weak genetic barriers to hybridization, so interspecific crosses can occur between genetically related taxa that co-occur and maintain overlapping flowering periods [7][8][9][10].Hybridization may influence species boundaries through differential introgression, local adaptation, speciation, or a partial merger of taxa termed reticulate evolution [3,10,11].Conversely, hybridization may negatively affect endangered species by disrupting co-adaptive gene clusters or eroding rare species [11,12].Orchid species delimitation is further complicated by symbioses [13,14].Biotic interactions with pollinators and mycorrhizal fungal distributions may affect the niche range of orchid taxa [15,16], form parapatric boundaries between sympatric orchid populations, or facilitate hybridization at sites of colocalization.Characterizing interspecific genetic exchange can uncover semipermeable species boundaries and help resolve morphology-based taxon delimitation.
Hybrid zones are the geographic area where genomic exchange occurs between diverged lineages in primary or secondary contact [17].Hybrid zone dynamics reflect the interaction of phenotypic divergence with adaptive selection, generating a plant evolution center that cannot be fully recapitulated under controlled settings [18][19][20][21].Conservation of orchid species in the light of hybridization depends on understanding the types of crosses between species, the intensity of genetic exchange, whether the hybridization dynamics lead to the exchange of advantageous or deleterious alleles, and whether hybridization is facilitated by anthropogenic activities [22].Intermediate phenotypes in overlapping orchid populations have historically been accepted as hybrids [23,24], but genetic data are needed to understand gene flow, intraspecific variation, and phenotypic plasticity [10,25,26].
The terrestrial orchid genus Platanthera (Rich.)[27] is a powerful system for studying plant hybridization because it maintains extensive diversity and species with overlapping spatial boundaries [19,25,28].Platanthera orchids span most of the temperate Northern Hemisphere and have undergone extensive floral radiation [10,29,30].Plastid and nuclear ITS markers showed evidence for polyploid hybridization and speciation in the Platanthera hyperborea complex [26].When neutral genetic markers (AFLPs) were used to genomically assess the sympatric zone between the closely related P. bifolia and P. chlorantha, researchers showed good separation between parent species but found intermediate plants maintained the genetic profile of P. bifolia [10].Finally, targeted enrichment of multilocus datasets reconstructed the phylogeny of Platanthera subgenus Limnorchis and was used to estimate that the subgenus diverged from Platanthera 3-4.5 Mya [30].With observed hybridization and backcrossing, Platanthera orchids may help evaluate the evolutionary significance of porous species boundaries.
Within Platanthera, species in subgenus Blephariglottis sect.Blephariglottis [31] demonstrate strong floral radiation [32,33].The Orange Fringed orchid, P. ciliaris (L.) Lindl, the Crested Orange Bog orchid, P. cristata (Michx.),and the White Fringed orchid P. blephariglottis (Willd.)Lind.have a wide distribution and overlap in several regions throughout the Eastern United States between New England and Florida and as far west as Texas [34,35].The three parent species are found in sphagnum bogs, pine savannas, meadows, and prairies, with the habitats of P. ciliaris being particularly varied [34].The specificity of two key biotic interactions, orchid mycorrhizae and pollinators, may constrain parent species distributions and generate parapatric boundaries between populations of species.While a range-wide analysis of pollinator fauna is missing for the subgenus Blephariglottis, studies have identified swallowtail butterflies (Papilio spp.) as pollinators of P. blephariglottis and P. ciliaris and bumblebees (Bombus spp.) as the predominant pollinator of P. cristata and an occasional pollinator of P. blephariglottis [34,[36][37][38][39][40].Partial overlap in pollinators and flower morphology may facilitate occasional cross-pollination events and interspecific gene flow [32].The distribution of specific mycorrhizal fungi partners impacts the distribution of parent species and whether hybrids can establish [15].Preliminary explorations of mycorrhizal diversity in P. ciliaris, P. cristata and P. blephariglottis have shown the three species to associate primarily with fungi in the genus Tulasnella [41][42][43], and an overlap in appropriate fungi might support the post-zygotic success of hybridization in sites of colocalization.In sites with species overlap, pairs of the three parent species have been observed to form putative hybrids with intermediate morphology [32].
Here we studied the putative crosses of P. blephariglottis with P. cristata to form P. x canbyi [33] and with P. ciliaris to form P. x bicolor [24] (Figure 1a).We analyzed sympatric sites where the hybrid's parent species co-occurred, sympatric sites where orange parent species co-occurred and P. blephariglottis was previously observed, and allopatric sites where a single parent species occurred (Figure 1b).Out of the ten sites studied, we found putative hybrids for P. x bicolor in Central Pennsylvania as well as coastal North Carolina, and P. x canbyi on the eastern shore of Maryland as well as New Jersey.The putative hybrids with intermediate phenotypes suggested parental species had crossed, but how they crossed and at what frequency was still unknown.The orange complex has been analyzed morphologically and observationally with pollinator ecotypes [32,34].Still, the assumed hybridization had not been confirmed with genomic data.An attempt to identify species and putative hybrids in the orange complex with standard DNA barcoding markers (ITS, trnH-psbA, trnL) failed to distinguish the taxa due to very limited sequence diversity.Sequences are available on NCBI GenBank (accession numbers OQ474558-OQ474566 and OQ507373-OQ507385).The genetic similarity across the hybrid complex for these markers likely reflects recent species divergence that is better resolved with multilocus analyses of genetic divergence [28,44,45].
To test the hypothesis that hybridization occurs between closely related Platanthera taxa, we conducted 3RAD genomic sequencing of allopatric and sympatric sites [44,45].We explored the extent of interspecific crosses of P. blephariglottis with P. ciliaris and P. cristata, and aimed to answer the following questions: (1) Is there hybridization and evidence of ongoing introgression?(2) What hybrid classes result from the interspecific genetic exchange?(3) How does genetic composition relate to morphological measurements?Are there cryptic hybrids?We then discuss how answers to these questions can inform strategic measures to preserve field sites with rare orchids and inform investigations into how anthropogenic change affects hybridization.

Hybridization and Backcrossing in the Fringed Orchids: Morphology
We compiled data on six morphological traits-lip width, lip length, lip with fringe, longest fringe, spur length, flower width, and inflorescence width-in the targeted populations from a larger dataset on variations in the complex [32] (Figure 1b) to assess morphological variation among taxa.We measured inflorescence width from the widest point of the inflorescence.We aimed for the center of the inflorescence for other morphological measurements, but we prioritized peak blooming flowers that occasionally fell at the bottom or top of the inflorescence, depending on the blooming stage.We obtained data on two flowers from each individual (n = 111) and according to flower color and site history, we assigned a putative species name of P. ciliaris (n = 23), P. x bicolor (n = 12), P. blephariglottis (n = 40), P. x canbyi (n = 10), or P. cristata (n = 26).We are missing flower width data for a P. blephariglottis sample in site PON (PON_ble_12) and spur length data for a P. cristata sample in SSS (SSS_cri_31).
We analyzed floral trait data in R 4.1.2[46].First, we calculated the mean and standard deviation of the traits for each putative species or hybrid with Morphotools [47].
To visualize the clustering of putative hybrid morphology, we implemented principal component analyses (PCA) of floral measures with MorphoTools and prcomp functions in R [47].We calculated permutational multivariate analysis of variance (Permanova) [48], a non-parametric approach that does not require prior information about the population distribution.We performed 9999 permutations of the permanovas and implemented Benjamini-Hochberg p-value correction for multiple comparisons.

Hybridization and Backcrossing in the Fringed Orchids: Genomics
For genomic analyses, we collected 147 samples from the Platanthera hybrid complex, composed of P. ciliaris (n = 32), P. x bicolor (n = 18), P. blephariglottis (n = 53), P. x canbyi (n = 9), and P. cristata (n = 37) in ten populations across four states between July and August of 2019 (Figures 1b and S8).We placed leaves in bags with silica gel immediately after collection for optimal preservation of DNA.We extracted DNA with the NucleoSpin Plant DNA extraction kit (Macherey-Nagel), following the manufacturer's instructions and using 20 mg of dried leaf material per sample and an elution volume of 50 µL.Three enzymes, MspI, BamHI, and ClaI, were used to digest the DNA, and RAD library preparations were conducted by the EHS lab at the University of Georgia [49].Paired-end 150 bp sequencing of pooled libraries was performed on an Illumina HiSeq instrument by Admera Health (South Plainfield, NJ, USA).
We conducted quality control validation on the raw sequence data with FastQC and removed twelve samples with fewer than 400,000 raw reads.We processed paired 3RAD reads from the remaining samples (n = 135) with iPyrad 0.9.81 [50,51] on the Smithsonian High-Performance Cluster (SI/HPC) [52].We demultiplexed each sequencing run with a maximum of one barcode mismatch and merged the demultiplexed reads for the remaining iPyrad steps.We increased the maximum depth to 1 ×10 6 within each sample to ensure coverage of clusters that were sequenced in high depth.We applied a clustering threshold of 90% minimum similarity because it is more appropriate for closely related species [53].We made the adapter filter stricter (2) to remove barcodes leftover from demultiplexing [51,54].We raised the maximum number of heterozygotes in the consensus to 0.1 to retain more alignments [53].Finally, we increased the minimum number of samples used per locus to 50 since downstream MrBayes and NewHybrids analyses are sensitive to missing data.The remaining iPyrad parameters were kept as default values [51].
To visually depict the degree of variation among the Platanthera taxa, we conducted dimensionality-reducing PCA.We filtered for SNPs that were present in 75% of samples (MinCov = 0.75), meaning SNPs not present in over 25% of samples were excluded [51].We then used the iPyrad analysis toolkit to perform PCA with a sample impute method and t-distributed Stochastic Neighbor Embedding (tSNE) plots of PC axes with a perplexity of 9 and 10 ×10 5 iterations [51].To model the evolutionary relationship between the Platanthera taxa, we conducted a Bayesian inference of phylogeny with MrBayes 3.2.7a on the Smithsonian HPC [51,55].Implementation of the Metropolis coupled Markov Chain Monte Carlo (MCMC) algorithm with four heated chains, 8 million generations, a burn-in fraction of 25%, and sampling of every 1000th tree constructed a phylogenetic tree that was visualized in Jupyter Notebooks with toytree [56].

Observed Hybridization Types
To visualize the allele frequencies of Platanthera taxa, we conducted admixture analyses with STRUCTURE [57].We ran ten runs each with K values 2-5 to conduct k-means clustering of the SNP datasets using a burn-in of 200,000 and 1 ×10 6 MCMC steps [58,59].We averaged values across runs with CLUMPP [60].To determine the best supported ancestral population value, we analyzed the highest rate of change in the log probability of data between successive K values [58,59].
We used NewHybrids to assign hybrid categories with posterior distributions from a Bayesian model-based clustering framework [61].First, we re-sorted the demultiplexed files from iPyrad step 1 into two sets.One set contained P. x bicolor and its respective parent species (P.blephariglottis and P. ciliaris).The second set contained P. x canbyi and its respective parent species (P.blephariglottis and P. cristata).We re-ran steps 2-7 of the iPyrad pipeline using the adjusted parameters explained previously.The P. ciliaris set (n = 92) contained 4443 loci and the P. cristata set (n = 84) contained 4927 loci.We implemented Fst analyses with hybriddetective to filter the variants into the 300 maximum allowed by NewHybrids to focus on the most informative loci between the two parent species [62,63].We converted the VCF file generated in iPyrad into a genepop file using PGDSpider with a population definition file that assigned the parent species to either Pop1 or Pop2 under the assumption that loci informative across parent species would be informative for the hybrids.We converted the Fst selected loci from a VCF file to a NewHybrids input file using PGDSpider [64] and generated hybrid assignments in parallelnewhybrids with a burn-in of 1 ×10 5 iterations, a sweep of 4 ×10 6 , and Jeffrey's theta and Pi priors [61,63,65,66].

Morphological versus Genetic Composition
We overlaid the genetic composition of the orchid populations on the morphological measurements we had visualized with PCA dimension reduction analyses.We used scatter-Pie [67] to replace the points in the morphology PCA with pie charts reflective of the genetic admixture of the population determined by STRUCTURE analyses.For samples paired across morphological and genomic analyses (n = 9), we populated a pie chart with the STRUCTURE output for the corresponding sample.For unpaired samples, we populated the pie charts with the averaged STRUCTURE output for the corresponding species or hybrid at the same sites.As we did not obtain genomic data for the P. cristata samples in site PON, we averaged STRUCTURE output of P. cristata in sympatric sites NAS and GSP.

Morphology Analysis
We found significant differences among species in flower morphology when we conducted a multivariate Permanova analysis of the combined traits (p < 0.001) and observed strong clustering in dimension reduction analyses (Figure S2) [32].The putative parent species, P. blephariglottis, P. ciliaris, and P. cristata, clustered independently.The measured floral traits from the putative hybrids, P. x bicolor and P. x canbyi, clustered between the appropriate parent species.Moreover, the two principal components described 73.8% and 17.8% of the variance in flower morphology (Figure S2).The sub-complex of P. blephariglottis, P. ciliaris, and its hybrid P. x bicolor maintained morphological similarity (Figures S1-S3).

Genomic Analysis
3RAD analysis yielded an average of 1.20 ×10 6 ± 6.9 ×10 5 reads per individual that maintained an average Phred quality score of 38.iPyrad analyses identified 149,669 single nucleotide polymorphisms (SNPs) across 7214 loci after filtering.PCA dimension reduction analyses of these SNP data sets showed clustering of the parent species and hybrids between the appropriate parent species, albeit with some overlap between P. cristata and P. x canbyi as well as the extensive overlap between P. ciliaris and P. x bicolor (Figures 2 and S4).The upper cluster of P. cristata overlapping with P. x canbyi were samples from sympatric sites where P. cristata was documented to have genetic admixture (NAS, NJT).The bottom P. cristata cluster contained samples from LI, GSP, and SSS.Similarly, the bottom cluster of P. ciliaris was composed of samples from the allopatric site WSP.Cluster patterns were maintained in tSNE plots (Figure S4).To investigate gene flow between the three parent species, we conducted posterior probabilities of the SNP datasets.An ancestral population size of K = 2 was the peak value in deltaK analyses [58], but the mean L(K) was highest for K = 3 and K = 4 (Figures S6 and S7).As the deltaK method is known for bias towards K = 2 [68] and K = 3 both had a high mean L(K) and made sense biologically, we presented plots for both K = 2 and K = 3.
The two clusters defined in K = 2 broadly correspond with a P. blephariglottis + P. ciliaris cluster and a P. cristata cluster (Figure 3), indicating gene flow is stronger between P. blephariglottis and P. ciliaris.Platanthera cristata was more distinct, with some signs of introgression in sympatric sites NJT and NAS and the hybrid P. x canbyi maintained varied admixture (Figure 3).K = 3 maintained the cluster corresponding with P. cristata and further divided P. blephariglottis and P. ciliaris into two clusters.These two clusters broadly correspond to current species perceptions of P. blephariglottis and P. ciliaris and reveal significant admixture between the two species, corroborating it as a K value worth exploring.
While the SNPs of P. cristata and P. blephariglottis in allopatric sites (BS, LI sites in Figure 3) were almost exclusively assigned to one ancestral population, allopatric P. ciliaris samples (WSP) maintained some gray-orange admixture (Figure 3).More admixture was observed in sympatric sites.Platanthera ciliaris in a sympatric site with P. blephariglottis (VB) and in two sympatric sites where P. blephariglottis was not found at the time of the study but was reported in the recent past (GSP, SSS) maintained considerable admixture.Sympatric P. x bicolor and P. ciliaris maintained nearly indistinguishable structure profiles.Admixture of putative Platanthera cristata and P. x canbyi samples varied across sites.Platanthera cristata samples in sympatric sites with P. blephariglottis (NAS,NJT) maintained a genetic admixture that was slightly stronger in site NAS relative to NJT.Moreover, P. x canbyi samples maintained split or red skewed admixture is sites NAS and PON.In sites where P. blephariglottis was only reported in the recent past (GSP, SSS), the P. cristata samples maintained negligible admixture comparable to the allopatric site LI.
We observed site dependency in the degree of genomic admixture between P. cristata and P. blephariglottis that reflected the current colocalization of parent species.The admixture profile of the hybrid complex reflected the geographic distribution of the parent species.Platanthera cristata is not found in Pennsylvania [35] and we found negligible red structure in the samples from the Pennsylvania site VB.The three parent species have been observed at a preserve in Brunswick, North Carolina.While extensive admixture was seen in P. ciliaris samples of this region (GSP, SSS), negligible admixture was seen in P. cristata (LI, GSP, SSS).Interestingly, some admixture was observed between P. ciliaris and P. cristata in the GSP+SSS region, indicating potential hybridization between the two orange species.Admixture across our sites reflected the absence of parent species in certain geographic regions, but reduced admixture was observed in most P. cristata and P. blephariglottis samples.Species were generally retained across the Bayesian inference trees (Figure 4).A clade with several subclades (A-D) contained exclusively P. blephariglottis samples.Another clade with two subclades (I,J) contained P. cristata and two paraphyletic groups of P. x canbyi (Ia,Ja).Several scattered clades contained P. ciliaris and the hybrid P. x bicolor (F,G,H), but allopatric P. ciliaris formed a distinct clade (E).We observed a site effect within the clades of parent species, as clade H only contained samples from sites in the same preserve (GSP and SSS), and the hybrids in NAS and VB clustered with their putative parents.We classified the posterior probability that an individual belonged to a certain genotype frequency class (F1, F2, backcross, pure parent).Samples from allopatric sites (BS, CRA, WSP, LI) were assigned to a pure parent class in NewHybrid analyses.Almost all of the putative P. blephariglottis samples were assigned to the pure parent genotype in both sympatric and allopatric sites.Most of the P. ciliaris samples in sympatric sites where P. blephariglottis was localized at the time of the study (VB) or recently localized (GSP, SSS) were identified as F2 P. x bicolor.The putative P. x bicolor hybrids were also predominantly F2 P. x bicolor, but a minority were identified as P. blephariglottis or backcrosses with P. blephariglottis (Figure 5).Platanthera cristata in sites where P. blephariglottis was recently localized (GSP, SSS) were all assigned parent genotypes.The putative P. cristata in sympatric sites where P. blephariglottis was localized at the time of the study (NAS, NJT) maintained a mixture of P. cristata parent species and backcrosses of P. x canbyi with P. cristata.Putative P. x canbyi samples also maintained a range in genotype assignment, with some found to be late generation F2 hybrids, backcrosses to P. cristata, and backcrosses to P. blephariglottis (Figure 5).We found that sympatric P. ciliaris and P. blephariglottis showed substantial signs of hybridization and that sympatric P. cristata and P. blephariglottis also hybridized, but to a lesser degree.Even in sympatric sites where P. blephariglottis was not found at the time of this study (GSP, SSS), most putative P. ciliaris were P. x bicolor hybrids, and none of the P. cristata were hybrids.To understand floral trait measurements with respect to genetic composition, we reconstructed the morphology PCA plot using pie charts composed of STRUCTURE admixture (Figure 6).We found that genetic-based name assignments maintained clear morphological trends, especially for P. blephariglottis and P. cristata.The samples on the outer fringes of the P. cristata and P. blephariglottis clusters were measured in sites where the genetic structure was relatively pure, and samples clustered in the center demonstrated more admixture.Platanthera ciliaris was more variable in genetic structure and morphology.The P. ciliaris cryptic hybrids maintained a range of morphotypes across the P. ciliaris distribution.Mor-phological or color-based techniques align in part with P. blephariglottis and P. cristata but do not model the genomic profile of P. ciliaris.

Hybridization in Allopatric and Sympatric Sites
Natural hybridization within the plant kingdom is a double-edged evolutionary force.Hybridization may facilitate adaptive allele sharing or break apart coadapted gene groups, just as it may generate stable hybrid populations or genetically erode rare parent species [12,69].We documented contemporary hybridization between the fringed orchid species P. ciliaris, P. blephariglottis, and P. cristata with genomic and morphological approaches.Our multilocus datasets documented interspecific gene flow across hybrid zones.We generated well-supported clades of P. cristata and P. blephariglottis but identified a convoluted phylogenetic distribution of P. ciliaris, P. x bicolor, and a subsection of P. blephariglottis samples.As we assigned species names based on color identification in the field, we did not expect clear species distinctions in the phylogenetic tree or a distinct clade of contemporary hybrids.Interestingly, we found Platanthera ciliaris and P. blephariglottis hybridized extensively, demonstrating weak reproductive barriers.Platanthera cristata hybridized to a lesser degree with P. blephariglottis in only sympatric sites.Although the three parent species are distributed throughout the Eastern United States, signs of hybridization generally occurred in regions with the appropriate parent species, suggesting that hybridization may be geographically limited by the distance between sites.
The limited rate of hybridization between P. blephariglottis and P. cristata might reflect additional site-specific prezygotic reproductive barriers, such as pollinator availability and flowering time, or post-zygotic reproductive barriers, such as seed inviability and mycorrhizal fungi constraining effects.Swallowtail butterflies (Papilio spp.) have been reported to visit P. ciliaris and P. blephariglottis and bees (Bombus spp.) have been reported to visit P. cristata and other species in the subgenus Blephariglottis [34,[37][38][39][40]).Some researchers argue that the smaller flowers and shorter spurs of P. cristata suggest that they are pollinated predominantly by bees [38,40], which may explain the more distinct morphology and genomic composition of P. cristata we observed.Differences in phenology also represent barriers to hybridization that could vary along the range of the species.Variations in flowering times across the distribution area of the species have not been assessed in detail, but during fieldwork, we observed a greater overlap in flowering time for P. cristata and P. blephariglottis in NJ and MD than in NC, agreeing with a higher rate of hybridization and admixture in the northern part of the sampling area for this species pair.Flowering time for P. ciliaris and P. blephariglottis overlapped in all sites where they co-occurred, consistent with more even rates of admixture in sympatric sites across the sampling area.The difference in morphology and especially spur length, between P. cristata and the larger flowered P. blephariglottis might also contribute to the lower rate of hybridization and admixture in this species pair.
If the two species overlap in mycorrhizal fungal partners, it is plausible that hybrids are compatible with strains used by both parents.The parent species have been found to associate with strains of fungi in the genus Tulasnella [41][42][43], and a single Tulasnella strain successfully germinated both P. blephariglottis and P. ciliaris [42].This suggests that mycorrhizal fungal specificity is likely not the primary postzygotic barrier to hybridization between P. ciliaris and P. blephariglottis, given that parental fungi are available.More research is needed to determine the mycorrhizal specificity of hybrids in relation to their parents.Whether differential OMF use by P. cristata and P. blephariglottis affects hybrid viability also remains to be tested.The strength of reproductive barriers between the two species varied but was fairly maintained across geographic areas, suggesting a balance between the reticulate gene flow within the hybrid zone and the preservation of parental genotypes [70].Selection against hybrids has been documented for Epidendrum neotropical orchids as they maintain prezygotic isolation and low hybrid fertility [6].Studies of several Orchis hybrid zones have also suggested that post-zygotic barriers are sufficient to maintain species integrity when reproductive isolation falters [70].Individual reproductive barriers are generally permeable, but species integrity is maintained by a combination of several reproductive barriers that arose through reinforcement [71,72], an important note for recently diverged species in this study.
Reproductive barriers and other ecological boundaries cap the extent of genomic exchange between diverged taxa.When minimal barriers are present, hybridization can disrupt genetic integrity or promote competition due to hybrid vigor [21,73].While P. blephariglottis and P. cristata maintained relatively limited hybridization, sympatric sites with P. ciliaris were primarily composed of P. x bicolor hybrids.We also found no P. ciliaris samples with negligible admixture in sympatric or allopatric sites, suggesting P. ciliaris might not be sufficiently isolated to prevent interspecific gene flow.Our sampling occurred over a large geographic area along the East coast of the United States, so we expect similarly high levels of hybridization to occur in other sympatric sites.While extensive hybridization documented between P. blephariglottis and P. ciliaris may support the two species being color morph varieties of one species [34], their genomic composition can largely distinguish the parent species and suggests they are closely related species that retain genetic mixing.Given our sampling was enriched in sympatric sites, future research into additional allopatric sites and multi-generational field studies can assess whether anthropogenic forces, such as agricultural expansion and erratic flowering, influence the observed hybridization and affect competition between P. x bicolor hybrids and their parent species [16,21,73].

Morphology versus Genetic Composition
The relationship between genomic and phenotypic variation has been tackled by key models, including Helianthus sunflowers, Heliconius butterflies, finches, and hamlet fish [63,74].The environmental factors governing the hybridization patterns of these models help define drivers of introgression and rapid speciation.Orchidaceae is among the most biodiverse flowering plant families.Still, studies characterizing the molecular diversity of orchids and their hybrids are limited compared to other plant families [75].As Platanthera orchids are a key angiosperm model for hybridization due to extensive adaptive radiation and selection mediated, in part, by specialized pollinator relationships [10,34], we coupled genomic and morphological approaches to characterize the genetic diversity of orchids.
Quantitative morphological characteristics partially reflect the hybrid complex, as putative hybrids with intermediate morphology tended to maintain mixed genetic admixture and hybrid assignment.While the extensive admixture and cryptic hybridization in P. ciliaris lacked any clear patterns with color-based or morphological measurements, the admixture and cryptic hybridization in the rest of our Platanthera complex corresponded with morphological measurements to a moderate degree.We speculate that the difficulty distinguishing the orange parent species and hybrids was partly because slight yellow is easier to distinguish from white than lighter orange is from solid orange, suggesting that color-based species identification of the white P. blephariglottis species is more precise than the orange species.
Field identification of the orange parent species (P.ciliaris and P. cristata) and their respective hybrids with P. blephariglottis based on flower color alone was not sufficient.As several visually identified P. cristata were identified as cryptic P. x canbyi, it is more appropriate to implement both color and quantitative morphological measurements when genomic analyses are impractical.These morphological measurements can focus on a small set of easily measurable traits, such as spur and fringe length [32], and be supplemented with genomic analyses when possible.The wide range of P. ciliaris morphological measurements were not reflective of genomic hybrid status, suggesting that genomic approaches are important to tease apart cryptic hybridization in the subcomplex P. ciliaris forms with P. blephariglottis.A limitation of our study was not pairing morphological and genomic measurements on an individual sample level.Direct sample pairing in future studies may allow better linking of morphology and genomics.The range of admixture profiles for the orange parent species demonstrates that the complex interplay between phenotypic divergence and molecular composition [18] in actively evolving hybrid zones is best resolved with genomic approaches.

Species Boundaries and Conservation Models
Hybridization and introgression influence orchid evolution [44], and the resulting genetic and morphological intermediates make it challenging to delimit species.The species boundaries of rapidly evolving Platanthera orchids undergoing interspecific gene flow demonstrate this complexity.A ploidy of 2n = 42 is broadly conserved across the majority of Platanthera taxa [76], and our results suggest repeated backcrossing supports P. x canbyi and P. x bicolor being homoploid hybrids and maintaining the same chromosome number as parent species [19].As the genomic structure of the Platanthera hybrid complex aligns with their modern species colocalization, the hybrids likely represent patchily distributed products of contemporary hybridization.Ongoing gene flow in hybrid zones is often a transitional phase within the longer-term interspecific exchange that can allow for variation in allele frequency [10,17].Unlike most hybrid zone models, however, Platanthera orchids do not occur along a single plane of colocalization.With patchily distributed species, the orange Platanthera complex maintained a continuum of genetic structure in sites with multiple parent species.Cross-pollination and the establishment of seedlings with the appropriate fungi is a differential process likely needed to generate stable hybrid populations [77] and may influence the patchy distribution of the Platanthera taxa.Interspecific exchange in the hybrid zones containing the Platanthera complex demonstrates continuity in species boundaries that merits future research into whether hybridization is driven by anthropogenic forces, including wetland drainage, agricultural expansion, invasive species expansion, pollinator decline as a result of habitat disruption and pesticides, and erratic flowering in a warming climate [16].
The three primary species models (biological, phylogenetic, and ecological) delimit species in different ways, but underlying evolutionary principles-such as variation, genetic determination, inheritance, selection, and accumulated change-are consistent.As speciation occurs as a continuum [3], it is arguably as important to analyze the evolutionary processes influencing diversification, such as hybridization.Natural hybridization may support biodiversity through novel breeding opportunities or impair hybrid fitness by generating genomic obstacles.Here, we outline a few potential impacts of our observed hybridization on species-level evolution.Hybridization localized to certain sites could support a local homogenization of parent species and potentially ephemeral speciation.
If the hybrids have an adaptive advantage over their parent species, they might evolve into a stable population with the potential to generate a separate species.Still, widespread hybridization can result in genetic erosion of endangered or locally endangered parent species [12,69].The hybrid P. x canbyi is confined to sympatric sites, and we do not see any clear evidence of genetic erosion in P. cristata or P. blephariglottis.However, the prominence of cryptic P. x bicolor might reflect the partial erosion of P. ciliaris.Conservation models may need to assess whether anthropogenic activities promote a breakdown of natural species barriers, especially between P. ciliaris and P. blephariglottis, to inform site preservation and whether to co-plant parent species in restoration efforts.If the processes leading to hybridization are natural, however, no intervention would be warranted.Longitudinal studies into whether anthropogenic factors, such as fragmentation, or natural factors, such as fungal specificity or species fitness, influence colocalization and promote species erosion would further inform conservation decisions.

Conclusions
To describe gene flow between possibly diverging taxa, we analyzed hybridization across a set of allopatric and sympatric sites.We showed that multilocus genomic analyses provided the resolution to study hybridization between closely related Platanthera taxa when single fragment approaches were insufficient.The genomic composition of the Platanthera hybrid complex revealed contemporary hybridization in sites of colocalization.We showed permeable reproductive barriers between P. ciliaris and P. blephariglottis and somewhat permeable reproductive barriers between P. cristata and P. blephariglottis.To determine whether anthropogenic activities drive breakdowns of natural barriers to gene flow, longitudinal studies into habitat fragmentation and symbiotic partner decline of the Platanthera hybrid complex may help explain the extensive hybridization we observed across the Eastern United States.Hybridization dynamics provide insight into the interplay between genomic structure and species boundaries.A genomic understanding of diversity and hybridization coupled with measures of anthropogenic pressure may inform Platanthera conservation.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/d15030384/s1, Figure S1 Funding: This research was funded by the Carlsberg Foundation (Grant: CF190076), the Garden Club of America (The Elizabeth Gardner Norweb Summer Environmental Studies Scholarship), the Washington Biologist Field Club, the Maryland Native Plant Society, the American Orchid Society, and the National Science Foundation (REU Award: #1950656).

Figure 1 .
Figure 1.(a) Platanthera hybrid complex.The parent species P. blephariglottis (gray) crosses with P. cristata (red) to form P. x canbyi (pink) and with P. ciliaris (orange) to form P. x bicolor (yellow).(b) Site acronyms represent the approximate locations of ten sites across the Eastern United States.Labels for sympatric sites (PON, NJT, NAS, GSP, and SSS) are colored black, and allopatric sites (WSP, CRA, LI, BS) are colored according to the taxon that was present.Circles on the map represent the geographic localization of the taxa.We assigned sites and taxa prior to genomic sequencing.

Figure 2 .
Figure 2. Principal Component Analysis of Platanthera taxa.PCA of SNP datasets reflected clustering of parent species based on their genomic composition.

Figure 3 .
Figure 3. Population structure of Platanthera taxa.Every SNP is assigned a probability of belonging to ancestral populations (K = 2 or K = 3) represented by gray, red, or orange.Assignments are summarized in a single bar.Below each bar is a symbol where color represents the field assigned species name, P. blephariglottis (gray), P. ciliaris (orange), P. cristata (red), P. x canbyi (pink), and P. x bicolor (yellow), and shape represents whether one (circle) or two (triangle) parent species were present at the site.

Figure 4 .
Figure 4. Phylogenetic relationship among Platanthera species and putative hybrids across sympatric and allopatric populations.The shapes correspond to sympatric sites where the hybrid's parent species co-occurred (triangle), sympatric sites where orange parent species co-occurred and P. blephariglottis was observed recently (square), and allopatric sites (circle).The Bayesian posterior probabilities of all deep branches were maximally supported (100), and the external branches also maintained high branch support (Figure S5).

Figure 5 .
Figure 5. Genotype assignment of the Platanthera hybrid complex.The bars are a cumulative probability for the different population types: pure parent 1 (gray), pure parent 2 (red for the top graph and orange for the bottom graph), F1 hybrid (pink for the top graph and yellow for the bottom), F2 hybrid (brown), backcross with parent 1 (blue), and backcross with parent 2 (purple).Below each bar is the sample identification with a color reflective of species.The allopatric sites with one parent species are BS, CRA, LI, and WSP.

Figure 6 .
Figure 6.PCA of floral morphology for genetically-identified Platanthera species.Each sample point shows a pie chart generated based on population structure data for the individual sample or an average of the same samples in that particular site due to discordance between the samples analyzed genetically and those analyzed morphometrically.For clarity, we circled the putative hybrids (pink and yellow) and the P. blephariglottis samples that clustered with P. ciliaris morphometrically (gray).