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

Abstract 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.


| INTRODUC TI ON
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 humanmodified 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 , increased competition from invasion (edge effects; , loss of symbionts (mycorrhizae, pollinators; Broadhurst & Young, 2006;Jacquemyn et al., 2004), reduced recruitment , 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(Aguilar et al., , 2008Angeloni 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 postfragmentation (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.

| Study species
The eastern prairie fringed orchid (Platanthera leucophaea) is a longlived, 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 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).
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,  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).

| 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 e ), 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 (F IS ) 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 eff , as it has good precision for microsatellite data with limited sample sizes and populations with small effective sizes (100-200; Waples & Do, 2008). N eff 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, 1997) and G ST which is equivalent to F ST 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)  with sample sizes larger than 25 using stratified. R code (https:// gist.github.com/mrdwa b/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 10 5 and collected data for 10 5 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 Rousset, 1997) and G ST (Pons & Petit, 1996) against Euclidian geographic distance for populations with sample sizes greater than 18 using the Mantel test in GENALEX (Peakall & Smouse, 2006

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 H e and N e as measures of genetic diversity, F is and relatedness as measures of inbreeding, and pairwise F ST and G ST , 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, N e , 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.
2.5.5 | 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 (N e ), expected heterozygosity (H e ), 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 (H e and N e ) and inbreeding (F is ) to the natural log of the average census size. This allowed us to look for trends that were less impacted by sampling bias. is a geographically disjunct population and likely an outlier, the next highest value was in Michigan (F is = 0.14).

| 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%;  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, t 33 = −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

| DISCUSS ION
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 easternto-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 F st 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) 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 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).

CO N FLI C T S O F I NTE R E S T
The authors do not have any conflicts of interest. Writing -review & editing (equal).