Anthropogenic fragmentation increases risk of genetic decline in the threatened orchid Platanthera leucophaea.

Protecting biodiversity requires an understanding of how anthropogenic changes impact the genetic processes associated with extinction risk. Studies of the genetic changes due to anthropogenic fragmentation have revealed conflicting results. This is likely due to the difficulty in isolating habitat loss and fragmentation, which can have opposing impacts on genetic parameters. The well-studied orchid, Platanthera leucophaea, provides a rich dataset to address this issue, allowing us to examine range-wide genetic changes. Midwestern and Northeastern United States. We sampled 35 populations of P. leucophaea that spanned the species' range and varied in patch composition, degree of patch isolation, and population size. From these populations we measured genetic parameters associated with increased extinction risk. Using this combined dataset, we modeled landscape variables and population metrics against genetic parameters to determine the best predictors of increased extinction risk. All genetic parameters were strongly associated with population size, while development and patch isolation showed an association with genetic diversity and genetic structure. Genetic diversity was lowest in populations with small census sizes, greater urbanization pressures (habitat loss), and small patch area. All populations showed moderate levels of inbreeding, regardless of size. Contrary to expectation, we found that critically small populations had negative inbreeding values, indicating non-random mating not typically observed in wild populations, which we attribute to selection for less inbred individuals. The once widespread orchid, Platanthera leucophaea, has suffered drastic declines and extant populations show changes in the genetic parameters associated with increased extinction risk, especially smaller populations. Due to the important correlation with risk and habitat loss, we advocate continued monitoring of population sizes by resource managers, while the critically small populations may need additional management to reverse genetic declines.

INTRODUCTION

Anthropogenic‐driven fragmentation is a major threat to biodiversity worldwide (Lienert, 2004). Over time, structural changes to habitat reduces patch size while increasing proximity to human‐modified landscapes and isolation between populations (Haddad et al., 2015). These interwoven landscape effects can operate over potentially long timescales to drive declines in both species and genetic diversity (Haddad et al., 2015; Ibáñez et al., 2014). Identifying the factors that most impact species decline can help managers prioritize conservation practices. There is debate over the relative importance of different types of habitat changes in driving biodiversity losses (Fahrig, 2017; Fahrig et al., 2019; Fletcher et al., 2018; Hadley & Betts, 2016). Although loss of habitat alone can explain biodiversity loss (Fahrig, 2017), patterns of biodiversity change are often explained by the complex interactions between patch area declines, connectivity reductions, and increased edge effects (Fletcher et al., 2018; Haddad et al., 2015). Given that genetic factors are one of the major causes of species extinction (Frankham, 2005), evaluating the impacts of anthropogenic changes on species success is critical to conservation efforts (Leimu et al., 2010).

The local extinction of a plant species is typically driven my multiple interacting factors. Under anthropogenic changes to landscapes such factors can include changes in habitat suitability (Breed et al., 2012), increased competition from invasion (edge effects; González‐Varo et al., 2012), loss of symbionts (mycorrhizae, pollinators; Broadhurst & Young, 2006; Jacquemyn et al., 2004), reduced recruitment (González‐Varo et al., 2012), loss of habitat area, and increased patch isolation (Butaye et al., 2001; Honnay et al., 2002). Together these changes negatively impact the demographic trajectory and genetic diversity, and ultimately reduce the reproductive output of plant populations (Aguilar et al., 2006, 2008; Angeloni et al., 2011; Honnay & Jacquemyn, 2007b; Leimu et al., 2006; Vranckx et al., 2011). These factors work together to accelerate decline, spiraling a population in a downward trajectory known as the “extinction vortex” (Gilpin & Soule, 1989). This increased risk is in part driven by genetic changes to the population during habitat fragmentation including increased inbreeding, loss of genetic variability, and increased divergence between populations (Lowe et al., 2005; Reed & Frankham, 2003). Together these genetic changes will have detrimental effects on population fitness and viability (Leimu et al., 2006) and influence the potential for a species to adapt to ongoing or future environmental changes (Jump et al., 2009; Lande & Shannon, 1996; Manel & Holderegger, 2013; Vilas et al., 2006).

The impact of fragmentation on a plant population can reduce pollen and seed dispersal, increasing genetic drift and inbreeding, while lowering diversity. Small populations have fewer potential mates, which can lead to increased bi‐parental inbreeding and even selfing. In time, this can lead to increased genetic load and inbreeding depression (Kramer et al., 2008). This increased genetic load can increase the risk of extinction if populations cannot withstand the consequences of inbreeding depression (Wallace, 2003). Increasing migration between small remnant populations, either through creation of corridors or movement of individuals, increases genetic diversity (Submitting author et al. in prep) and is used as a management tool for populations impacted by fragmentation (Pavlova et al., 2017; Whiteley et al., 2014). However, in some plant species, there is evidence that fragmentation can have a positive impact on gene flow (Breed et al., 2015; Matesanz et al., 2017), although it has been suggested that this paradox is the product of sampling individuals established pre‐fragmentation (Breed et al., 2013; Vranckx et al., 2011), rather than seedlings or younger individuals that arose post‐fragmentation (Breed, Marklund, et al., 2012; Breed et al., 2013; Breed, Stead, et al., 2012; Vranckx et al., 2011; Yates et al., 2007).

Successfully predicting the effects of anthropogenic habitat loss and fragmentation on genetic diversity and species resilience requires an integrated approach that considers genetics, demographics, isolation, and habitat quality across a species’ range (Lowe et al., 2005). To this end, we conducted a range‐wide genetic analysis that included assessment of population size, isolation, and habitat availability of the federally threatened orchid Platanthera leucophaea (Nuttall) Lindley (eastern prairie fringed orchid). We sampled populations spanning from critically small (n < 5) to large, robust population sizes (n > 1000). Our objectives were to (i) investigate range‐wide genetic patterns in P. leucophaea to determine if we can detect the underlying genetic structure of the species before fragmentation, (ii) determine the relationship between population size and genetic variation, and (iii) determine if changes in genetic parameters might indicate if populations impacted by fragmentation are heading toward extinction.

METHODS

Study species

The eastern prairie fringed orchid (Platanthera leucophaea) is a long‐lived, perennial, terrestrial orchid (Bowles, 1983, Figure 1) native to high‐quality wetland communities including wet prairies, sedge meadows, fens, and bogs of North America. The species was once found in contiguous wetland communities east of the Mississippi River, south of Ontario, and north of Kentucky. Historically, populations numbered in the thousands (Bowles, 1983; COSEWIC, 2003; Figure 2); however, due to habitat loss and fragmentation, the total number of populations has decreased by over 70% in the United States (USFWS, 1999). Remaining populations occur in small remnants with diminished numbers of individuals (from 1 to over 1000 individuals; Table1). As a consequence, P. leucophaea was federally listed as threatened in the United States in 1989 and federally listed as endangered in Canada in 2005 (COSEWIC, 2003). Today, there are 97 documented populations of P. leucophaea known to exist in the United States (Illinois, Indiana, Iowa, Maine, Michigan, Missouri, Ohio, Wisconsin (USFWS, 2016)) and 21 in Canada (Ontario, COSEWIC, 2003).

FIGURE 1

Platanthera leucophaea, the eastern prairie fringed orchid. Photo: Rachel Wells

FIGURE 2

Present day and historic range map of P. leucophaea in North America, with populations sampled in blue, all other extant populations in black, and locations of herbarium specimens in orange

TABLE 1

Genetic data for 26 large and 10 critically small populations (n < 18) of Platanthera leucophaea that span the range of the species. Population metrics include name code (abbr.), state, latitude and longitude, sample size, and census size in 2015. The genetic measures calculated include effective population size (Neff), average number of alleles (Na), effective number of alleles (Ne), genetic diversity (He), inbreeding coefficient (Fis), relatedness (r), average pairwise Fst, and average pairwise Gst. Contact the author for location information

Code	State	Samples size	Census size (2015)	Eff. Pop. Size (Neff)	# Alleles (Na)	# Eff. Alleles (Ne)	Genetic div. (He)	Inbreeding (Fis)	Relatedness (r)	Ave. Pairwise Fst	Ave. Pairwise Gst
MS	IA	30	75	38.40	7.63	4.82	0.74	0.12	0.04	0.08	0.08
BM	IA	30	600	56.80	9.00	6.15	0.79	0.11	0.04	0.09	0.09
CB	WI	29	75	58.30	7.88	5.03	0.67	0.07	0.03	0.08	0.08
FM	WI	30	63	44.80	7.63	4.39	0.74	0.11	0.04	0.07	0.07
KW	WI	30	100	62.60	9.00	6.11	0.77	0.15	0.02	0.06	0.07
UPA	WI	30	100	33.80	8.00	5.09	0.76	0.05	0.04	0.07	0.07
CH	WI	31	75	95.60	8.88	4.81	0.71	0.14	0.02	0.06	0.06
BS	IL	29	30	112.90	7.38	3.71	0.70	0.08	0.08	0.09	0.09
WT	IL	30	33	38.40	9.00	5.36	0.71	0.18	0.02	0.06	0.06
MC	IL	30	39	341.40	7.88	5.31	0.75	0.08	0.04	0.08	0.07
LY	IL	30	103	57.60	8.25	5.52	0.75	0.07	0.03	0.06	0.06
GC	IL	30	110	12.50	6.63	4.57	0.70	0.12	0.08	0.10	0.10
W	IL	30	347	43.70	9.13	5.35	0.72	0.13	0.02	0.07	0.07
SOM	IL	30	394	27.10	6.88	4.03	0.69	0.20	0.07	0.08	0.08
PM	MI	18	18	434.20	7.38	4.47	0.68	0.04	0.04	0.08	0.07
SB	MI	30	100	112.80	7.75	4.12	0.66	0.10	0.06	0.10	0.09
LCL	MI	31	38	40.50	6.75	3.43	0.67	0.07	0.13	0.15	0.13
MO	MO	30	57	107.10	7.38	4.08	0.69	0.08	0.08	0.10	0.09
YA	OH	23	23	183.60	7.38	4.36	0.70	0.01	0.06	0.08	0.08
MTZ	OH	25	25	82.20	7.75	4.88	0.74	0.07	0.06	0.09	0.08
PC	OH	28	334	Inf	8.38	5.34	0.72	0.06	0.03	0.06	0.06
NT	OH	29	135	46.50	8.50	5.58	0.76	0.08	0.03	0.06	0.06
YO	OH	29	1316	36.90	6.38	3.90	0.70	0.12	0.09	0.12	0.11
KB.83	OH	30	68	115.40	8.13	4.23	0.73	0.01	0.04	0.09	0.08
LD	OH	30	132	Inf	8.88	5.41	0.74	0.19	0.01	0.06	0.06
CC	OH	30	654	81.60	8.00	4.92	0.74	0.15	0.05	0.07	0.07
Critically small populations
Long	IL	6	6	12.70	4.50	3.72	0.69	−0.34	0.10	0.10	0.08
IN	IN	5	5	9.70	4.50	3.68	0.68	−0.09	0.01	0.07	0.06
ME	ME	5	15	3.70	2.63	1.87	0.39	0.44	0.20	0.32	0.28
BG	MI	10	10	78.50	5.63	3.96	0.66	0.01	0.08	0.08	0.08
DC	MI	7	7	31.20	4.63	3.33	0.59	−0.04	0.13	0.22	0.18
SW	MI	13	158	25.60	5.38	3.58	0.66	0.14	0.07	0.10	0.10
MI	MI	7	7	11.70	4.00	2.83	0.60	−0.22	0.15	0.15	0.13
MDW	OH	13	15	22.50	6.88	4.49	0.72	0.09	0.03	0.07	0.07
KB.N	OH	6	8	63.80	3.38	2.49	0.52	−0.05	0.17	0.17	0.15
KB.S	OH	3	6	1.80	3.25	2.85	0.63	−0.47	0.09	0.07	0.06

Platanthera leucophaea can self‐pollinate but relies on pollinators to produce seed. The main pollinators are nocturnal flying hawkmoths in the family Sphingidae (Bowles, 1983; Cuthrell et al., 1999). P. leucophaea flowers produce dust‐like seed that is thought to be wind dispersed. Hence, despite the highly fragmented distribution of extant populations, long‐distance pollination by hawkmoths and seed dispersal by wind may facilitate some gene flow between populations. Germination requires mycorrhizae (Zettler & Piskin, 2011) before a protocorm is formed, which may stay dormant for several years depending on nutrients supplied by mycorrhizal fungi. Plants may take 2–13 years to reach reproductive maturity (USFWS, 2016).

Study sites

For this study, we selected 36 P. leucophaea populations that span its US range (Figure 2) and vary in population size, patch isolation, composition, and area. In states that the species is extant, we targeted populations that were historically the largest (based on available census sizes from the USFWS; Table S1). From each population, leaf tissue was haphazardly collected from 30 flowering individuals or all flowering individuals, if a population had less than 30 plants. One population (ME) had no flowering plants and only vegetative plants were sampled from this population. We collected 4–5 cm from the leaf tip of each plant, and dried the tissue in silica gel for later DNA extraction. We sampled 25 populations in 2015 and 3 in 2016. DNA extractions provided by Lisa Wallace (Wallace, 2002) were used to fill in sampling gaps in Michigan (one population, 1999) and Ohio (three populations in 1998). We were particularly interested in including those populations that are critically small (<50) and therefore most likely to be impacted by inbreeding and loss of diversity associated with an extinction vortex. One of the challenges of including these populations is that the sample sizes are below the number recommended for accurate genetic assessment (Hale et al., 2012). To address this, we generated a sampling effort curve to determine the minimum sample size needed to give equivalent results; similar to species accumulation curves used in ecology (Fisher et al., 1943).

Molecular data

To characterize the genetic structure of each population, genomic DNA was extracted using a CTAB extraction protocol (Doyle & Doyle, 1987). DNA quality was estimated using a NanoDrop 2000 spectrophotometer (Fisher Scientific). Thirty‐one microsatellite primers previously developed by Ross et al. (2013) for Platanthera praeclara were screened for genotyping of P. leucophaea. Of the 31 primers, 9 did not amplify in P. leucophaea, 10 produced monomorphic peaks and 12 produced consistent banding patterns. Twelve microsatellite loci were identified for genotyping P. leucophaea, which amplified and produced consistent and reliable banding patterns (PP05, PP07, PP09, PP10, PP16, PP19, PP22, PP23, PP24, PP27, PP30, and PP31; Ross et al., 2013). These were used for genotyping all individuals using fluorescently tagged forward primers (Sigma‐Proligo). PCR reactions were performed in a 10 µl reaction mixture containing the following: 3 µl DNA, 0.5 µl forward and reverse primers, 0.125 µl BSA and 0.875 µl H20, and 5 µl PCR master mix 2x (Promega, Madison). PCR was run at 94°C for 5 min, then 35 cycles of 94°C for 40s, 60°C for 40s, 72°C for 3 min, and then a final extension at 72°C for 10 min for six primers (PP05, PP07, PP22, PP16, PP19, and PP27), while for the remaining primers (PP09, PP10, PP23, PP24, PP30, and PP31), the extension step was extended to 1 min. Genotypes were scored using a CEQ 8000 Genetic Analysis System and CEQ FRAGMENT ANALYSIS software (Beckman Coulter).

Population variables

Population size, patch composition, and degree of geographic isolation were measured or calculated for all populations in multiple ways. Using multiple measurements is particularly important for population size, as orchid numbers can fluctuate significantly from year to year depending on environmental and local biotic conditions (USFWS, 1999). As metrics for population size, we included: (1) a census of population size at time of leaf collection in 2015 and 2016; (2) a calculated median population size, as well as minimum and maximum size in populations from yearly censuses (2003–2015); and (3) the categorical census size determined by the USFWS in the recovery plan (0 ≤ 10, 1 = 10–25, 2 = 25–50, and 3 > 50; USFWS, 1999; Table S1).

The metrics for patch composition were calculated using a 1, 10, and 20km buffer around each population and tabulating the landcover‐type areas within, using classifications of the 2011 Gap Analysis Program (GAP), National Landcover Data. As the 1, 10, and 20km proportions were all equivalent, we ultimately focused on the 1 km buffered area. To identify the spectrum of preferred habitat types of this species, we generated a list of the natural landcover types which intersected each population centroid from across the sampling area. We assumed, given that the species was found within these landcover types, that they represent potentially suitable habitat (Table S2). Therefore, to calculate the total potential patch area within the 1km buffer around each population, we summed the area for all landcover types in which P. leucophaea was found. The remaining landcover types were then characterized as either natural, agricultural, or developed landcover types (Table S2). Natural area was defined as any landcover which was not developed, agricultural, or suitable habitat for P. leucophaea. Once all landcover types were classified, we summed the area to give us the total area of suitable habitat, natural area, agricultural, and development within the 1km buffer surrounding each population. Ultimately, the metrics associated with patch composition included (1) the USFWS ranked categorical habitat sizes (0 = <2.5 acres; 1 = 2.5<62.5 acres; 2 = 62.5<125 acres, and 3 = >125 acres; USFWS, 1999), (2) the amount of suitable habitat (patch size), (3) natural area, (4) agricultural land, and (5) total development.

The degree to which each population was isolated was determined using the average pairwise distance between all populations, average distance to the 5 and 10 nearest extant populations, and landscape resistance. The average Euclidian pairwise distance was calculated in SPAGeDi (Hardey & Vekemans, 2002) from the latitude and longitude of each population. The average nearest‐neighbor distances to the 5 and 10 closest neighbors were calculated in ArcGIS 10.3.1 using all known extant populations. Landscape resistance is a distance metric calculated using Circuitscape, which accounts for the variability in landscape types for movement between populations. This required categorizing all of the landcover types from the Multi‐Resolution Land Characteristics consortium, National Land Cover Database (NLCD), from 0 to 5 based on its suitability to pollinator movement or habitat suitability for colonization (Table S3). Circuitscape then tests multiple trajectories through the landscape and produces a relative metric of the distance between populations based on landscape suitability and determine the shortest distance with the least landscape resistance connecting two populations (Shah & McRae, 2008).

Statistical analysis

Genetic parameters

We used 12 microsatellites to measure factors commonly impacted by fragmentation and associated with declining populations which included: (1) genetic diversity (Schlaepfer et al., 2018), (2) effective population size (England et al., 2010), (3) inbreeding levels within (Leimu et al., 2010; Schlaepfer et al., 2018), and (4) differentiation between populations (Miles et al., 2019). All primers were tested for departure from Hardy–Weinberg Equilibrium (HWE) at the locus, population, and global levels, using Genepop (Raymond & Rousset, 1995). The potential of null alleles and mis‐scoring was tested using exact tests in Micro‐Checker (Van Oosterhout et al., 2004). Genetic diversity was quantified using effective number of alleles per locus (N), and expected heterozygosity (H e), which were calculated in GENALEX (Peakall & Smouse, 2006). Both of these metrics are less sensitive to differences in sample size than other metrics. To measure rates of inbreeding within a population, we calculated Weir and Cockerham’s (1984) estimates of Wright’s inbreeding coefficient (FIS) and Queller and Goodnight’s Relatedness (R) in GENALEX (Peakall & Smouse, 2006). We used the linkage‐disequilibrium method in NeEstimator V2.01 (Do et al., 2014) to measure the population size, N, as it has good precision for microsatellite data with limited sample sizes and populations with small effective sizes (100–200; Waples & Do, 2008). N is an important parameter because it is not always reflected in inbreeding and genetic diversity variables (Bazin et al., 2006; Raymond & Rousset, 1995) and is impacted by fragmentation (England et al., 2010). Finally, to measure genetic differentiation, we used Rousset’s linearized FST (FST/(1e−FST; Rousset, 1997) and GST which is equivalent to FST but more appropriate for microsatellites (Pons & Petit, 1996). The pairwise genetic and spatial distances were calculated in SPAGeDi (Hardey & Vekemans, 2002).

Sample size

Since the sample sizes in 12 populations were below the accepted sampling size (n < 25) for attaining accurate population genetic metrics (Hale et al., 2012), we generated a sample effort curve for He and Fis, similar to that proposed by Bashalkhanov et al. (2009) to identify the minimum sample size which produces accurate estimates of genetic parameters. To achieve this, we subsampled all populations with sample sizes larger than 25 using stratified. R code (https://gist.github.com/mrdwab/6424112). We randomly sampled individuals from each population to generate populations that ranged in sample sizes from 2 individuals to 25 individuals and repeated 10 times for each sample size. We then calculated the expected heterozygosity using the R package StrataG.R (Archer et al., 2017) and calculated inbreeding using the R package adegenet (Jombart, 2008) for all subsamples. To determine the minimum sample size at which heterozygosity and inbreeding results plateaued, we plotted the mean standard error of all samples using ggplot2 (Wickham 2016; Supplemental 3).

Range‐wide genetic structure

To determine if there were range‐wide patterns in genetic structure within P. leucophaea, we used the Bayesian clustering analysis software STRUCTURE v2.2 (Falush et al., 2007; Pritchard et al., 2000) to visualize population subdivision (number of genetic clusters, K) among our 36 study populations. This software uses individual multilocus genotypes to test for the presence of population structure without a priori assignment of individuals to populations by finding population groupings with the least possible disequilibrium using a Markov Chain Monte Carlo method. We carried out 40 independent runs per K using a burn‐in period of 105 and collected data for 105 iterations for K = 1–40. The minimum value of K that can explain the data was assessed using the rate of change in the log‐likelihood probability of data between corresponding K values (ΔK) as detailed in Evanno et al. (2005) using Structure Harvester (Earl & vonHoldt, 2012).

To identify isolation by distance, we compared pairwise genetic distances FST (FST/(1−FST); Rousset, 1997) and GST (Pons & Petit, 1996) against Euclidian geographic distance for populations with sample sizes greater than 18 using the Mantel test in GENALEX (Peakall & Smouse, 2006). To determine if other genetic variables varied by location, we tested for correlations of HE and NE (genetic diversity), FIS and relatedness (inbreeding), and Neff (effective population size) against latitude and longitude using the correlation coefficient panel for pairs function in R Statistical Software (R Core Team, 2016).

Modeling populations and genetic parameters

Principal component analyses

We used the princcomp function to create a principal component analysis (PCA) summarizing all variables into a single metric that captures the spectrum of variation for that trait. All variables were first tested for normality and were log‐transformed when appropriate (including area of habitat, developed area, natural landcover, census size in 2015, and the minimum, maximum, and median population sizes across the available years of the USFWS censuses (IL and OH only)). For population size, we used five variables: census in 2015 (Census), median census size (MedCensus), minimum and maximum size (Min & MaxCensus), and the categorical census size (CatPopSize). For patch composition, we used four variables: percent suitable habitat, percent natural area, percent agriculture, and percent development. Finally, for patch isolation, we used four measurements: average pairwise distance (Distance), average distance to 5 and 10 nearest extant populations (NeN5, NeN10), and landscape resistance as calculated in Circuitscape (resistance). For each PCA, we used the get_pca_var function in the R package factoextra 1.0.6 (Kassambara & Mundt, 2020) to determine the amount of variation explained by each axis (Eigenvalue and proportion of variation explained) and the contribution of each metric in explaining the spread. The predict function in Vegan (Okansen et al., 2019) was used to extract the coordinates for each population along the PCA axis. To check for independence of our different model parameters, we used cor.test function in R to look for association between our PCA axes and latitude and longitude.

Models

To investigate whether population and landscape variables explain genetic parameters, we used linear models that included all of the extracted axes that cumulatively explained at least 80% of the variation for the population size, patch isolation, and patch composition PCAs (see Table 2 for list of each metric used), as well as latitude and longitude to account for geographic variation in these traits. The simplest model was selected through backward and forward elimination using the StepAIC function in R Statistical Software. The best model was compared for homoscedasticity and then tested against the null hypothesis using the ANOVA function. The genetic parameters tested in the model were He and Ne as measures of genetic diversity, Fis and relatedness as measures of inbreeding, and pairwise FST and GST, which are measures of fixation index often used as a proxy for genetic distance between two populations. As some populations were too small to accurately reflect genetic metrics of diversity and inbreeding, we created two datasets, one that contained data from all populations and a second that was restricted to just larger populations (n > 18). In addition to genetic parameters, we were also interested in how population sizes vary with anthropogenic changes. Therefore, we modeled both census size and effective population size, Ne, against the PCA axes calculated to estimate the impacts of anthropogenic change, patch isolation, and patch composition. Since effective population size is not always correlated with census size, we modeled this term with and without the size PCA axes. After running the models, we found that the two disjunct populations in Missouri and Maine were two‐ to fourfold further from any other populations, making them large geographic outliers that had a disproportionate impact on the model analyses. Since these were effectively isolated, and to get a more accurate representation of processes in the center of its range, we excluded these populations from the model analyses.

TABLE 2

Breakdown of PCA results for a) population size, b) patch area, and c) patch isolation, including eigenvalues, proportion of variance explained, and breakdown of percent explained by each variable used in the analysis

PCA	All data (36 populations)				Only larger sample size (n > 18; 26 populations)
a) Population Size	pcaSize1	pcaSize2	pcaSize3	pcaSize4	pcaSize1	pcaSize2	pcaSize3	pcaSize4
Eigenvalue	3.23	1.11	0.32	0.26	2.36	1.50	0.59	0.47
Total Variance Explained	65%	22%	7%	5%	47%	30%	12%	10%
Breakdown by Parameter
Census	24.05	1.55	15.03	57.46	24.59	2.80	31.75	39.09
MinCensus	6.97	68.17	0.54	0.03	6.75	54.31	1.03	0.10
MaxCensus	19.94	23.17	9.72	18.43	16.28	30.77	2.64	25.98
MedCensus	26.30	5.44	1.36	22.56	29.48	10.51	0.06	24.31
CatPopSize	22.74	1.67	73.35	1.53	22.90	1.61	64.52	10.53

Impact of small population sizes on genetic data predictions

To determine if the trends in genetic parameters associated with including the critically small populations is a product of sampling bias or decreasing population sizes, we randomly selected seven individuals (the average size of our critically small populations) from all populations sampled using the stratified.R package. From this subset of data, we recalculated the effective number of alleles per locus (Ne), expected heterozygosity (He), and Weir and Cockerham’s (1984) estimates of inbreeding in GENALEX (Peakall & Smouse, 2006). We repeated this 10 times to generate multiple measurements of each parameter. These 10 replicates of each population were combined into a single dataset. We then used linear models in R to compare genetic diversity values (He and Ne) and inbreeding (Fis) to the natural log of the average census size. This allowed us to look for trends that were less impacted by sampling bias.

RESULTS

Descriptive statistics of populations

Four of the 12 microsatellite primers were excluded (PP09, PP10, PP16, and PP19) from our analyses after testing for null alleles, departure from Hardy–Weinberg equilibrium (HWE), and mis‐scoring. The sample effort curve revealed that the minimum sample size needed for an accurate measurement of genetic diversity (He) and inbreeding (Fis) for our dataset was n = 18. For those populations (>18 individuals), the average number of alleles by population (Na) ranged from 6.3 to 9.1 (ave 7.9), the number of effective alleles (Ne) ranged from 3.4 to 6.1 (ave 4.8), and expected heterozygosity (He) ranged from 0.66 to 0.79 (ave of 0.66). For the 10 populations where sample sizes were below 18, the range of the number of alleles (Na) was from 2.6 to 6.8 (ave 4.5), the number of effective alleles (Ne) ranged from 1.8 to 4.5 (ave 3.3), and expected heterozygosity (He) ranged from 0.39 to 0.72 (ave of 0.61), all of which were at the lower end of the ranges seen. The range for Weir and Cockerham’s (1984) estimates of Wright’s Fis varied from 0.01 to 0.20 (Fis = 0.10), suggesting a moderate level of inbreeding across populations; however, populations less than 18 spanned a large extreme from highly outcrossed (Fis = −0.47) and highly inbred (Fis = 0.44). Although, the highest measured inbred population was located in Maine (Fis = 0.44), which is a geographically disjunct population and likely an outlier, the next highest value was in Michigan (Fis = 0.14).

Analysis of population genetic structure

Structure Harvester identified two genetic clusters (K) between our populations with a slight geographic gradient from west to east (Supplemental 1). Interestingly, the geographically disjunct population in Maine was comprised of almost equal amounts of both clusters. Pairwise genetic relatedness slightly increased with distance; the pairwise genetic distances (Fst) for population pairs greater than 18 individuals ranged from low (0.01) to relatively high (0.16) with a low average of 0.06 (Supplemental 2). Nei’s D, pairwise genetic distance, followed a similar trend, increasing slightly with distance between populations.

Principal component analyses

Differences in patch composition of all populations was 100% explained by the first four PCA axes (Table 2). The first axis (pcaArea1), which explained 43% of the variation, described the degree of development within the patch, being positively associated with the log area of urban development (41%) and negatively associated with the log area of agricultural development (38%). The second axis (pcaArea2), which explained 28% of the total variation, showed a strong association with total area of preferred habitat types (80%; Table 2; Figure 3). The final two axes explained 19% and 10% of the remaining variance, respectively (Table 2). For populations with >18 individuals (large populations), the first two axes explained most of the variation (71%; Table 2).

FIGURE 3

PCA variables and their contributions to the (a) habitat PCA, (b) isolation PCA, and (c) size PCA for all Platanthera leucophaea populations sampled

Differences in patch isolation between only the large populations was 100% explained by the first four axes, although the first axis explained 89% of that variation. Each metric of fragmentation (landscape resistance, average pairwise distance to 5 nearest and 10 nearest extant populations, and geographic distance) contributed equally to the variation along the first axis (Table 2). PCA axis 1 explained 81% of the differences in patch isolation between all populations with equal contributions from all four metrics (Table 2; Figure 3).

To quantify population size, we used four variables (census in 2015, median census size, minimum and maximum census size, and the categorical population size). Differences between population sizes were 98% explained by the first four PCA axes using four measurements of population sizes (census in 2015, median census size, minimum and maximum census size, and categorical population size), with 47% of the variation explained by census 2015, median census size, and categorical population size along axis 1 (pcaSize1). The second axis explained an additional 30% of the variation and was explained by the minimum census size. The remaining two axes explained 12% and 10% of the variation (Table 2). For all populations, the first axis explained 65% of the variation and was equally influenced by all metrics except minimum census size. By contrast, the second axis explained an additional 22%, which was mainly explained by minimum census size (68%; Table 2; Figure 3).

Correlations between model parameters

Before we ran the models, we tested all PCA axes for evidence of correlation between axes, as well as between latitude and longitude. Most pairs of axes showed no significant correlations with each other. An exception was patch isolation (pcaIso1), which was significantly negatively correlated with development (pcaArea1), and this was regardless of whether we used the dataset that includes all populations (r = −0.40, t 33 = −2.5, p = .02) or only those with larger sample sizes (r = −0.40, t 25 = 2.18, p = .04). And, patch isolation (pcaIso1) was also significantly negatively correlated with patch area (pcaArea2) regardless of whether we used the dataset that include all populations (r = −0.52, t33 = −3.5, p = .001) or only those with larger sample sizes (r = −0.57, t 25 = 3.45, p = .02). This suggests that for our dataset we were unable to parse out completely the impacts of patch size (pcaArea2) and development (pcaArea1) from the degree of patch isolation. Interestingly, in this data, increasing isolation was associated with decreased development and larger population sizes. To investigate if there was a geographic pattern to the fragmentation parameters, we tested for correlations with longitude and latitude. We found no correlations between longitude and size (pcaSize1 or pcaSize2), but there was a significant positive relationship with patch size (pcaArea2) for both the complete dataset with all populations (pcaArea2, r = 0.58, t = 4.13, p = .0002) and the dataset restricted to only larger populations (pcaArea2, r = 0.62, t = 3.95, p = .0005), suggesting that patch areas to the east are larger. Similarly, the significant negative correlations between longitude and patch isolation (all pops, r = 0.34, t = −2.13, p = .04; and large only, r = 0.44, t = −2.48, p = .02) suggest most populations to the eastern edge of the range are also less isolated. Additionally, there was a significant negative relationship between latitude and patch isolation (all pops, r = −0.47, t = −3.10, p = .004; and large only, r = −0.50 t = −2.93, p = .007), suggesting southern populations are less isolated. There were no correlations between latitude and patch size (pcaArea2), development (pcaArea1) or size (pcaSize1), but there was a significant positive correlation with the second PCA axis of size (pcaSize2) and latitude for both datasets (all pops, r = −0.38, t = −2.05, p = .05; and large pops only, r = 0.42, t = 2.66, p = .01). Although this axis only explained 6–12% of the variation for size, it is strongly negatively influenced by minimum census size. This might suggest that northern populations have recorded some of the smallest census sizes.

The relationship between genetic parameters and population metrics

For the models that included all populations, the variables that best explained the effective number of alleles (Ne; R 2 = 0.53, R 2 adj = 0.49, F 3,30 = 11.55, p < .001) included population size (pcaSize 1, p < .001 & pcaSize 2, p = .15) and longitude (p = .03 & plarge only = 0.17; Table 3). Similarly, the best model for expected heterozygosity (R 2 = 0.60, R 2 adj = 0.53, F 5,29 = 8.5 p < .001) also retained size (pcaSize1, p < .001) and longitude (p = .01), but also included patch isolation (pcaIso1, p = .13), amount of development (pcaArea1, p = .19), and latitude (p = .04). When we used the more restricted dataset the best model for effective number of alleles (Ne) was less supported (R 2 = 0.23, R 2 adj = 0.16, F 2,23 = 3.5, p = .05) but still included both population size (pcaSize 1, p = .07) and longitude (p = .17). While with expected heterozygosity, the best model (R 2 = 0.46, R 2 adj = 0.36, F 3,22 = 4.59, p = .01) also retained size (pcaSize1, p = .16) and longitude (p = .005), but it also included latitude (p = .01). In all scenarios, the genetic diversity metrics decreased with population size (Figure 4) and surprisingly longitude (from east to west), which is likely associated with greater development in the northwest of the range. The models that best explained levels of inbreeding (Fis) in the full dataset (R 2 = 0.41, R 2 adj = 0.40, F 1,32 = 22.5 p < .001) and the large population‐only dataset (R 2 = 0.38, R 2 adj = 0.26, F 4,21 = 3.2, p = .03) both retained population size a predictor (pcaSize1, p < .001 and p = .14, respectively), and the large populations dataset also retained development (pcaArea1 p = .04) and longitude (Longitude p = .02) as strongest predictors. By contrast, the models that tested for Relatedness (R) using only large populations did not retain any variables, but when we included all populations, the best model (R 2 = 0.32, R 2 adj = 0.29, F 2,31 = 7.62, p = .001) was negatively correlated with population size (pcaSize1, p = .001) and patch area (pcaArea2, p = .05). The models for degree of genetic differentiation (Fst, Gst) when using only large populations did not retain any variables; however, when we included all populations, the models for both Fst (R 2 = 0.28, R 2 adj = 0.21, F 3,30 = 3.97, p = .02) and Gst (R 2 = 0.26, R 2 adj = 0.18, F 3,30 = 3.55, p = .03) retained population size (pcaSize1, p = .009 and p = .02, respectively), patch isolation (pcaIso1, p = .09 and p = .08, respectively), and patch area (pcaArea2, p = .04 and p = .04, respectively). In all cases, the degree of differentiation increased with decreasing population size (Figure 4), decreasing patch area, and increased patch isolation (Table 3).

TABLE 3

Results of the linear models (F statistic (F), degrees of freedom (df), p value (p) and regression fit (both multiple and adjusted R squared)) comparing population characteristics, including population size (pcaSize1 and pcaSize2), patch isolation (pcaIso1) and patch area (pcaArea1, pcaArea2, and pcaArea3), and latitude and longitude against the following genetic parameters: effective population size (Neff), genetic diversity (Na, Ne, and He), inbreeding (Fis and relatedness), and differentiation (Fst and Gst). We also tested how population size is dependent on patch isolation, patch area, latitude, and longitude. Models were repeated with (a) all data (36 populations), and (b) restricted to only larger (sample size (n) > 18) sample sizes (26 populations)

Variable	F	df	p	Multiple R 2/ AdjustedR 2	pcaSize1	pcaSize2	pcaIso1	pcaArea1	pcaArea2	Lat	Long
(a) All data (36 populations)
Na	29.67	1,32	<.001	0.48/0.47	<0.001	‐	‐	‐	‐	‐	‐
Ne	11.55	3,30	<.001	0.53/0.49	<0.001	0.15	‐	‐	‐		0.03
He	8.53	5,28	<.001	0.60/0.53	<0.001	‐	0.13	0.19	‐	0.04	0.003
Fis	22.55	1,32	<.001	0.41/0.40	<0.001	‐	‐	‐	‐	‐
Selfing	NS	‐	‐	‐	‐	‐	‐	‐	‐	‐	‐
Relatedness	7.62	2,31	.002	0.32/0.29	0.001	‐	‐	‐	0.05	‐	‐
Fst	3.97	3,30	.02	0.28/0.21	0.009	‐	0.09	‐	0.04
Gst	3.55	3,30	.03	0.26/0.18	0.02	‐	0.08	‐	0.04	‐	‐
Census Size	3.04	2,31	.07	0.16/0.10	NA	NA	‐	‐	0.02	‐	0.07
Effective Population Size	3.35	1,32	.08	0.09/0.06	NA	NA	‐	0.08	‐	‐	‐
(b) Large sample size (n > 18; 26 populations)
Na	NS	‐	‐	‐	‐	‐	‐	‐	‐	‐	‐
Ne	3.50	2,23	.05	0.23/0.16	0.07	‐	‐	‐	‐	‐	0.17
He	4.59	4,21	.01	0.46/0.36	0.16	‐	‐	0.12	‐	0.01	0.005
Fis	3.2	4,21	.03	0.38/0.26	0.14	‐		0.04	‐	0.11	0.02
Selfing	4.57	3,22	.01	0.38/0.30	0.01	‐	0.01	‐	‐	0.07	‐
Relatedness	NS	‐	‐	‐	‐	‐	‐	‐	‐	‐	‐
Fst	NS	‐	‐	‐	‐	‐	‐	‐	‐	‐	‐
Gst	NS	‐	‐	‐	‐	‐	‐	‐	‐	‐	‐
Census Size	1.99	1,24	.17	0.07/0.03	NA	NA	‐	‐	0.17	‐	‐
Effective Population Size	3.29	2,23	.05	0.22/0.15	NA	NA	‐	‐	0.16	‐	0.01

FIGURE 4

Expected heterozygosity (He) (a–c), inbreeding (FIS) (d–f), and differentiation (FST) (g– i) by PCA size axis 1 (a, d & g), PCA isolation axis 1 (b, e, & h), and PCA patch area axis 2 (c, f, & i) for all populations of Platanthera leucophaea. Large populations, n > 18, are displayed with black circles, and small populations are hollow circles. Linear model (all populations), (a): p ≤ .001; (b): p = .13; (c): p = NS; (d): p ≤ 0.001; (e): p = NS; (f): p = NS; (g): p = .009; (h): p = .09; and (i): p = .04). Solid line, p ≤ 0.05, dashed line, p > .05, p = NS, not supported in linear model

The models for census size using only large populations (R 2 = 0.07, R 2 adj = 0.03 F 1,24 = 3.29, p = .05) retained only patch area (pcaArea2, p = .17), while the model using all population (R 2 = 0.16, R 2 adj = 0.10, F 2,31 = 3.0, p = .07) retained both patch area (pcaArea2, p = .02) and longitude (p = .07). For effective population, only the models using large population (R 2 = 0.22, R 2 adj = 0.15 F 2,23 = 3.28 p = .05) retained patch area (pcaArea2, p = .16), while models using all populations (R 2 = 0.09, R 2 adj = 0.07, F 1,32 = 3.4 p = .07) retained development (pcaArea1, p = .08) and longitude (p = .01). This suggests that patch area is the best predictor of population size, with larger patches having larger populations. The association of size and longitude likely reflects the higher number of small populations to northeast of the range (Table 3).

Impact of small sample size on models

For critically small populations, sample size can make it difficult to draw conclusions when comparing genetic data from populations with larger samples sizes (>30). To correct for differences in sample sizes, we subsampled all populations to generate populations with seven randomly selected individuals (hereafter, all populations). We found that the adjusted average values of NE, HE, and FIS measured on only seven individuals decreased compared to the original value (Supplemental 4). The magnitude of change somewhat explains the difference in the range of values seen between the 27 larger population (>18) and 10 smaller populations (≤18). For genetic diversity metrics, NE and HE, the average standard deviation across the 10 replicates was within the range of the difference seen between subsampled and complete datasets (st.dev = 0.61 and 0.06, respectively). By contrast, the average standard deviation between the complete and subsampled data for inbreeding (FIS) was much greater (SD = 0.15). The large range of variation in the inbreeding coefficient likely reflected that populations comprised of only seven individuals will vary more in the ratio of related and unrelated individuals. Although we found negative FIS values within subsampled datasets for most populations, all replicates of the critically small populations’ subsampled datasets had low‐to‐negative inbreeding values suggesting that for these populations this result is not a product of chance or small sample size alone.

Despite a clear downward trend in absolute values in all three metrics for the subsampled datasets, we found that there was a significant negative relationship between the log of census size and both measures of genetic diversity (NE, F 1,368 = 113.6, p < .001; and HE, F 1,368 = 98.36, p < .001). Interestingly, FIS of the subsampled populations was also negatively influenced by population size based on the 2015 census size (F 1,368 = 43.32, p < .001), like the results from the complete population samples (Figure S4). Hence, although the absolute values from the subsampled data were lower, the negative trend persisted (Figure S4). We interpret this to mean that the values for these smaller populations are a product of both sampling size bias and biological changes associated with small population sizes.

DISCUSSION

In this range‐wide study of the federally listed eastern prairie fringed orchid, Platanthera leucophaea, we found low levels of differentiation between populations, and moderate levels of inbreeding and genetic diversity within populations. There was a weak positive correlation between genetic and geographic distance from the eastern‐to‐western end of the orchid range (approx. 2000 km). Southern populations had slightly more genetic diversity than the northern populations. The low genetic differentiation and weak isolation by distance suggest historically high levels of gene flow between P. leucophaea populations. Over the last 100 years, human development and agricultural encroachment have fragmented P. leucophaea habitat, and as a consequence many populations have gone extinct and over 70% of the remaining populations have been reduced in size (USFWS, 1999). Of the many changes due to fragmentation, population size was consistently the best predictor of all genetic factors. For other parameters, the degree of genetic differentiation was shown to increase with patch area and isolation, while increasing patch isolation and development were associated with lower diversity and increased inbreeding. These relationships were particularly evident for the critically small populations (<18 individuals remaining). Therefore, a good indicator of potential genetic issues, when fragmentation occurs, is a drop in population size associated with less habitat area, especially when in an urban matrix (Miles et al., 2019). Given the yearly fluctuations in population sizes in orchid species, consistent monitoring is an important method to assess for extinction risk (Mace & Purvis, 2008).

Across the extant range of the species, we found relatively high levels of genetic diversity in many populations. The high diversity seen across the range is comparable to high levels of genetic variation found by previous studies (Havens & Alex Buerkle, 1999; Wallace, 2002). Many of the metrics used for diversity showed a slight trend with latitude and longitude and suggest that there is higher diversity to the southwest and lower diversity to the northeast. Considering the likely post‐glacial migration routes of P. leucophaea, the higher genetic diversity in the southwest could also be a product of it representing historic refugia (Hampe & Petit, 2005), which is supported by the lower diversity in the northeast of its range (Paul et al., 2013). A recent study in a sister species, P. praeclara (Ross & Travers, 2016) observed a similar range of diversity to P. leucophaea but higher number of alleles per locus, which might be a product of ascertainment bias, since these markers were developed in P. praeclara.

Despite the fragmentation, we found that larger populations show low levels of differentiation and low to moderate inbreeding. This is not surprising given that this orchid was historically widespread, occurred in large populations (USFWS, 1999), and is pollinated by long‐distance fliers (Bawa, 1990; Haber & Frankie, 1989; Linhart & Mendenhall, 1977; Skogen et al., 2019). These results differ somewhat from previous genetic studies, but this difference is likely driven by their focus on smaller populations, which also show elevated differentiation values compared to the larger populations in our study (Havens & Alex Buerkle, 1999; Paul et al., 2013; Wallace, 2002). Studies in the closely related P. praeclara, with similar hawkmoth pollinated flowers and wind‐dispersed seed, found comparable levels of genetic differentiation to this study (Pleasants & Klier, 1995; Ross & Travers, 2016). Our findings are also consistent with average levels of genetic differentiation found in other species of Orchidaceae (Machon et al., 2003; Phillips et al., 2012). Although P. leucophaea is now rare and populations are highly fragmented, the observed low Fst across the range is consistent with what is expected for a once historically widespread perennial species with wind‐dispersed seeds and long‐distance pollinators (Duminil et al., 2007; Loveless & Hamrick, 1984).

Population size consistently had the strongest relationship with changes in genetic diversity and inbreeding within a population. Past studies in P. leucophaea did not find that population size, measured as either harmonic mean (Wallace, 2002) or census size of the year sampled (Havens & Alex Buerkle, 1999), was correlated with genetic variation. This is likely an issue of sampling a limited number of populations and the high variability that most measurements for population size will produce. Orchid populations can fluctuate dramatically from year to year (Shefferson et al., 2019), making it difficult to produce a measurement that accurately reflects the population status. In our study, we overcame these issues by increasing the number of populations and using a PCA to incorporate different measures of population size (census in 2015, median, and the categorical population size) and fluctuations (census minimum and maximum size). Given that smaller populations are more susceptible to genetic drift and elevated inbreeding, the patterns we saw follow expectation, especially when considering the critically small populations. Previous literature reviews have also found a similar relationship to the size and population genetic diversity regardless if the species are rare or not (Honnay & Jacquemyn, 2007a). More importantly, P. leucophaea population size was positively correlated with a loss in fitness as well as genetic diversity (Leimu et al., 2006), supporting the concerns that these factors interact to accelerate the extinction vortex.

Surprisingly, we found that as P. leucophaea population size decreased, inbreeding also decreased, contrary to previous reviews which found a negative or no correlation (Angeloni et al., 2011). The positive relationship in our populations is likely driven by two biological phenomena: 1) that large populations of P. leucophaea maintain moderate‐to‐low levels of inbreeding, and 2) most of the populations with critically low numbers had negative inbreeding coefficients. Previous studies have reported similarly high levels of inbreeding for larger populations of P. leucophaea (Paul et al., 2013; Wallace, 2002) and P. praeclara (Ross & Travers, 2016). High inbreeding values in orchids are not unexpected. Orchids package their pollen into pollinaria, which ensures greater transport efficiency and maximizes the number of pollen grains reaching the stigma to assure high seed set (Harder, 2000; S. D. Johnson et al., 2005). However, the percentage of pollinaria that reach the stigma of another flower is low (<5% in Disa (S. D. Johnson et al., 2005)); roughly half of individuals rely on self‐pollination (Nilsson et al., 2009; Luyt & Johnson, 2001; Tremblay, 1994; Peakall, 1989; Salguero‐faria & Ackerman, 1999; Nilsson, 1992; S. D. Johnson et al., 2005). The other reason for this positive relationship is the negative inbreeding values in critically small populations, which suggests that mating is disassortative (i.e., only between unrelated individuals; Rasmussen, 1979; Stoeckel et al., 2006). Several different mechanisms have been proposed to explain this phenomenon, especially in critically small populations. The first is mate limitation, where species that cannot self will favor selection for the rarer allele (or sex; Hoebee et al., 2007; Sujii et al., 2015). Similarly, in species that are clonal (Halkett et al., 2005) or capable of unisexual reproduction (Johnson & Jonathan Shaw, 2015), there can be strong linkage disequilibrium generating large negative inbreeding values as the genetic lines diverge. Given that P. leucophaea is self‐compatible and non‐clonal, we do not think this is the situation in our study. The other possibility is that selection is favoring heterozygosity (heterosis). Although heterosis is beneficial in some systems (Stilwell et al., 2003; Oakley et al., 2015; Oakley et al., 2019; Bensch et al., 2006; W. H. Lowe et al., 2017), the fact that we are only seeing this in the smallest populations suggests that selection against homozygosity, associated with elevated inbreeding depression, is more likely (Aguilar et al., 2019; Angeloni et al., 2011; Charlesworth & Charlesworth, 1987).

The reduction in fitness associated with elevated homozygosity, known as inbreeding depression, is commonly seen in orchids (Juillet et al., 2007; Ortiz‐Barney & Ackerman, 1999; Sletvold et al., 2012). Several studies have shown a reduced seed production in selfed plants, including orchids (Sletvold et al., 2012), the genus Platanthera (Gregg, 1990; Nilsson, 1983; Patt et al., 1989; Travers et al., 2018), and more specifically this species (Wallace, 2003). However, inbreeding depression can also express itself in other demographic life stages including germination, seedling performance, and survival (Juillet et al., 2007; Sletvold et al., 2012). It is not uncommon for inbreeding coefficients to decrease as seeds transition to seedlings and then into adults, when most genetic sampling occurs (Aguilar et al., 2019; Cabin et al., 1998; Del Castillo, 1994; Honnay et al., 2008; Oostermeijer et al., 1995; Richards, 2000; Tonsor et al., 1993). The one critically small population that had elevated inbreeding levels was Maine, but it is also the only population where all sampled individuals were juvenile, supporting the hypothesis of higher inbreeding values in earlier life stages (Tonsor et al., 1993). Thus, we hypothesize that the negative inbreeding coefficient observed in our smaller populations is a result of a limited number of individuals surviving to flowering due to loss of individuals to inbreeding depression in early stages. This has been observed in other orchid species (Juárez et al., 2011) and will likely threaten the future viability of these populations (Spielman et al., 2004), emphasizing the need for conservation efforts within small P. leucophaea populations.

Additionally, we found support that important predictors of change to inbreeding, genetic differentiation, and diversity were patch isolation and the degree of urban development. The strength of these relationships was weaker than the population size, although this is likely because the three predictors were highly correlated in opposite directions. Many of the smaller populations were surrounded by greater urbanization but were also closer together and therefore less isolated. Thus, it was harder to untangle the impact of patch isolation versus development and patch area on our populations. This was reflected in our model of census size, which was associated with longitude and patch area, where the smaller populations on the western edge of the range had less available natural habitat. Similarly, the effective population size was lowest where urban development was the highest. Populations which showed some of the largest drops in population size over the censused years (minimum census size) were those populations to the north. Together these analyses show that populations in the Chicago Metropolitan area, at the northwestern edge of the range, have some of the smallest population sizes. The association between urbanization and population size suggests these populations are at high risk of going extinct.

The lack of relationship between genetic parameters and patch isolation suggests that it has less of an impact in this species. However, the lack of association in our data may be driven by many of our small populations occurring in close proximity, hence preventing us from disentangling the impacts of population size from patch isolation. However, as we also saw a lack of isolation by distance and weak differentiation across the range, this is likely a legacy of high rates of gene flow, at least in the past. This lack of differentiation is surprising, due to a long history of fragmentation and short generation times, but has been seen in other long‐distance pollinated species (Breed et al., 2013; Skogen et al., 2019). Long‐distance gene flow events, although commonly thought to be rare and difficult to document, are important for maintaining diversity in many systems. Hawkmoths in particular are of interest due to the fact that unlike bats and birds, which travel large distances but typically return to a nesting or roosting site daily, they are long‐distance flyers that migrate and are not home‐site specific. For these reasons, hawkmoths are expected to contribute to long‐distance gene flow despite multiple factors that may impede gene flow of P. leucophaea populations in modern fragmented landscapes (Aguilar et al., 2019).

Platanthera leucophaea has life history traits that are often associated with successful range edge expansion into newly hospitable habitat (long‐distance dispersal, highly mobile pollinators, and self‐compatibility; Parmesan, 2006). However, the loss of available habitat and reduction in population sizes will likely have long‐term consequences for the species. We found that a number of populations are showing signs of genetic decline, with evidence that populations are suffering from inbreeding depression and loss of genetic diversity. This was most evident in urban areas where the smaller populations have restricted patch area. Anthropogenic changes in patch area are directly related to population size changes, and we found that population size was a good indicator of genetic changes to populations of this orchid. Therefore, monitoring of these populations should continue to be prioritized in order to avoid population extinction due to genetic decline. More specifically, we note that populations with less than 15 flowering individuals are of highest concern and may lead to a demographic bottleneck if left unmanaged. Successive years of low census size may create a smaller “effective population size” and can therefore have serious genetic consequences. For some populations, which are showing signs of genetic decline, the augmentation of these populations, either through seed or pollen addition, is likely warranted. As many of these small populations have small neighboring populations, it is possible to conduct genetic augmentation with low risk of outbreeding depression, especially given the low differentiation recorded in this species (Amos et al., 2002; Frankham et al., 2011; Ralls et al., 2018).

CONFLICTS OF INTEREST

The authors do not have any conflicts of interest.

AUTHOR CONTRIBUTION

Claire Ellwanger: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Laura Steger: Investigation (equal); Writing – review & editing (equal). Cathy Pollack: Conceptualization (equal); Writing – review & editing (equal). Rachel Wells: Investigation (equal); Writing – review & editing (equal). Jeremie Fant: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal).

Supplementary Material

Click here for additional data file.