Distribution, phenology, and land cover preference risk assessments of four South-Central US milkweeds (Asclepias spp.) important for monarch butterflies


 South-Central US milkweeds (Asclepias spp.) are critical adult nectar and larval food resources for producing the first spring and last summer/fall generations of declining eastern migratory monarch butterflies (Danaus plexippus). MaxEnt niche models were developed for North American ranges of four important South-Central US milkweeds: Asclepias asperula ssp. capricornu, A. viridis, A. oenotheroides, and A. latifolia. Twelve models per species utilized subsets of six to eight of 95 edapho-topo-climatic variables chosen by a random subset feature selection algorithm. Milkweed weekly phenology was compared between early and late season periods of monarch activity. Novel land cover preference risk assessments were developed for milkweeds through land cover utilization-availability analyses, incorporating a novel sample bias reduction method for citizen science data before calculation of relativized electivity index (E*) land cover preference. Asclepias a. ssp. capricornu and A. viridis occurred more frequently during early season monarch activity, while A. oenotheroides and A. latifolia occurred more frequently during late season monarch activity. Milkweed utilization of roadsides varied from 6–31%. Developed-Open Space and Grassland Herbaceous land classes generally had highest benefit among milkweeds. Cultivated Crops and Shrub/Scrub had high risk. Combined milkweed high Ei* kernel density estimation surfaces resolved interior and coastal corridors of milkweed land cover preference providing functional connectivity for the monarch spring and fall migrations. A potentially critical gap in milkweed land cover benefit connectivity was identified in South Texas. Milkweed land cover preference risk assessments can be used to prioritize milkweed habitat conservation for enhancing monarch migration connectivity across the South-Central US.

use of glyphosate herbicides among corn and soybean crops (Pleasants 2017;Stenoien et al. 2018). Waterbury et al. (2019) found that A. speciosa, A. incarnata, and other Northwestern US milkweeds rarely utilized land covers of cultivated crops, developed, pasture/hay, and mixed and evergreen forest, and major habitat threats included invasive plants, herbicide use, and mowing. In comparison to data on A. syriaca in the Midwest (Hartzler 2010;Pleasants 2017;Kaul and Wilsey 2019) and A. speciosa in the Northwest, information is lacking on the utilization and risk of different LULC types for South-Central US milkweeds. We take a hierarchical approach (see Wilting et al. 2010 and Online Resource 1 Introduction) to assess species land cover preference risk for milkweeds, rst de ning the core habitat of milkweeds using an ensemble of feature-selected MaxEnt niche models. We then use a subset of high spatial resolution milkweed occurrences from citizen science iNaturalist data within the core habitat to develop novel land cover preference risk assessments for milkweed habitat from land cover utilization/availability analyses. Ogawa and Mototani (1991) performed the only other study of which we are aware that applied utilization/availability analysis to assess plant land cover selection (for further background on utilization/availability analyses, see Online Resource 1 Introduction). We also develop a novel citizen science data sample bias reduction method to reduce the weight of milkweed observations in oversampled areas of higher human population density and closer road proximity. Land cover classi cations for milkweeds were assessed using both a land cover GIS database and manual classi cation with publicly available internet aerial and ground imagery, and milkweed percent utilizations of roadside right of way (ROW) versus non-ROW habitats were also determined.
The objectives of this study were to (1) develop MaxEnt niche models that identify core edapho-topo-climatic habitats for each of the four study species, A. a. ssp. capricornu, A. viridis, A. oenotheroides, and A. latifolia, (2) compare the relative weekly phenology of the four milkweed species with respect to early and late season monarch activity within modeled South-Central US milkweed core habitats, and (3) use land cover utilization/availability analysis and a land cover electivity index to develop novel individual and seasonal combined land cover preference risk assessments across the South-Central US. Results are discussed in the context of monarch conservation.

Materials And Methods
Occurrence data Historical (1834Historical ( -2018 milkweed occurrence data throughout the entire North American ranges for the four study species were obtained from iNaturalist (2018a), SEINet (2017) David Berman and Kristen Baum, Oklahoma State University [OSU], 2017). Additional data representative of the known county distributions of milkweed species (Kartesz 2015) was obtained through georeferencing SEINet, vPlants, and OVPD records in various counties from which we otherwise lacked speci c locations. From 388 to 979 occurrence points per species were available for modeling after spatial thinning to 10 km. Partitioned data of subspecies for wide ranging species can yield better quality niche models than data of the combined subspecies (Gonzalez et al. 2011). Consequently, we partitioned data for A. a. ssp. capricornu to the East from A. a. ssp. asperula to the West (for details on the subspecies boundaries, see Online Resource 1 Methods). For milkweed phenology analyses, 2011-2021 iNaturalist milkweeds records over MaxEnt core habitats (see below) in the South-Central US were thinned to remove exact date/year/location duplicates, and further thinned to within 30 km for every year/week combination to reduce spatiotemporal bias in phenology. In utilization/availability analyses for the land cover preference risk assessment, 2014-2018 (within two years of 2016 land cover data) iNaturalist milkweed occurrence data at ca. 30-m resolution within core habitat for the South-Central US study area were spatially thinned to 1-km (for all occurrence data Excel les, see Online Resource 1 Table 1).

Study area boundaries
The background evaluation extents for niche modeling each of the four focal South-Central US milkweed species were represented by 500-km buffers around convex hull polygons from the milkweed occurrence points, which together extended over much of central eastern, southeastern, and southwestern North America. The land cover preference risk assessments for each milkweed species were con ned to a South-Central US study area comprising Arkansas and the majority of Oklahoma, Texas, and Louisiana ( Fig. 1). This study area excluded Mexico to the South due to the lack of a land cover data set for Mexico that matched the thematic resolution of the 2016 National Land Cover Dataset (NLCD) (Homer et al. 2015, Multi-Resolution Land Characteristics Consortium [MRLC] 2017) for the US, especially in terms of distinguishing important Pasture/Hay and Grassland/Herbaceous land cover classes (Fig. 1). The western edge of the South-Central US region consisted of the western distribution limit for A. a. ssp. capricornu in West Texas and Oklahoma (Woodson 1954) (see Occurrence data, above). The eastern South-Central US boundary consisted of the eastern boundary of Arkansas and northern Louisiana, and the Mississippi River in southern Louisiana (Fig. 1).

Niche modeling
An initial set of 95 environmental variables at 1-km resolution was used in developing MaxEnt niche models, including 57 climatic indices (Hijmans et al. 2005;Zomer et al. 2007Zomer et al. , 2008Trabucco and Zomer 2010;WorldClim 2017), 14 topographic indices (Lehner et al. 2008;Evans et al. 2014), and 24 edaphic indices (Hengl et al. 2017) (Online Resource 1 Table 2). Anthropogenic and land cover variables were intentionally excluded from the niche models to facilitate a more independent assessment of species a liation with land cover through the land cover preference risk assessment (see below) (for justi cation of initial variable selection, see Online Resource 1 Methods). MaxEnt (Phillips et al. 2006) version 3.3.3 models were developed from milkweed locations using the Rsoftware (R Core Team 2017) dismo package (Hijmans et al. 2011). The 1-km milkweed resolution occurrence locations were thinned by a 10 km spatial lter for reducing sample bias and spatial autocorrelation among the data (Boria et al. 2014). We generally followed methods described in Tracy et al. (2018) for generation of background and pseudoabsence points and calculation of AUC, TSS, and over tting (for details, see Online Resource 1 Methods).
For each milkweed species, we utilized the random subset feature selection algorithm (RSFSA) to evaluate performance of thousands of MaxEnt models developed from random combinations of subsets of 3 to 25 of the 95 variables, limiting correlation of variables to less than |0.7| (Tracy et al. 2018 Phenology relative to monarch activity Spatiotemporally thinned iNaturalist milkweed occurrence data with plant images were used to develop weekly stacked phenophase histograms and phenoperiod Gantt charts for niche modeled core habitat of each milkweed in the South-Central US. Phenophase histograms were visually assessed to estimate the best weeks for distinguishing between early and late season growth, which for all species was regarded as the periods from 1-33 weeks (January to mid-August) and 33-52 weeks (mid-August to December). The early and late season plant frequency data were tted to curves using CurveExpert software (Hyams 2010). Milkweed species weekly phenology was analyzed with respect to two important 11-week periods of monarch activity in the South-Central US (unpublished data) based upon monarch larval observations from Journey North ( Land cover preference risk assessment The hazard that various 2016 NLCD land covers represented for milkweed colonization were categorized using the estimated preference of each milkweed species for each land cover class within the South-Central US MaxEnt identi ed core habitat (see above)The spatial combination of hazard (milkweed preference) and exposure (milkweed core habitat) was used to develop milkweed land cover preference risk assessment maps. Thus, our risk assessment represents land use preferences of milkweed species and not risk as measured by direct land cover-related hazards or threats (e.g., herbicide use; . Milkweed land cover preferences were estimated based upon a utilization/availability analysis and an electivity index for occurrence of milkweeds among 2016 NLCD land cover classes. Land cover utilization/availability analysis relies heavily on accurate land cover classi cations for milkweed occurrences. Images for individual iNaturalist milkweed occurrences were examined to determine the NLCD land cover class within a 15 m radius of the location and within the context of the adjacent surroundings according to the 2016 NLCD legend (MLRC 2021) using historical imagery for 2014-2018 from Google Earth Pro (2021) and Google Maps (2021) Street View. Manually determined land cover classi cations for milkweed occurrences were compared with classi cations derived from the NLCD 2016 spatial database. In addition, the general habitat of the milkweed location was manually classi ed using the imagery into road ROW versus non-road habitats, and habitat utilization was analyzed among species using Chi-squared tests for independence (for further details on manual land cover classi cations using imagery and further general habitat subdivisions, see Online Resource 1 Methods).
The observed proportions of milkweeds in various land cover classes are also sensitive to errors from sample bias in citizen science data related to human population density and distance to roadways (e.g., Geldmann et al. 2016, Zhang andZhu 2018). Milkweed occurrence counts were weighted to reduce these biases using a novel target group-based Histogram Bin Ratio Weighting (HBRWt) method. This method assumes that (1) the spatial distribution of a large target group of taxa related to the taxa of interest represents a good estimate of citizen science sample bias from various sources (e.g., Ranc et al. 2017), and (2) the percent frequency distribution of the target group occurrences with respect to the spatial sources of bias (e.g., road proximity), should ideally approximate the percent frequency distribution of randomly generated points with respect to the sources of bias. These approximating assumptions provide an objective basis of correcting for sample bias that leads to more conservative estimates of milkweeds occurring within human population-dense areas or in proximity to roads. The HBRWt method involved rst constructing histograms of either human population density or distance to roads for two sets of data in the study area: (1) 10,000 random points; and (2) up to 10,000 target group points, where the target group is 30m resolution 2014-2018 iNaturalist records for major families of monarch preferred nectar plants (e.g., Asteraceae and Apocynaceae). In a particular human population bin range, if the random point percent frequency was 10% and the target group percent frequency was lower at 5%, the HBRWt adjustment factor for population density, PopDenHBRWt, was 10%/5% = 200% upwards for any milkweed point occurring within the population density bin (a raw count of one became two). The road distance HBRWt correction factor, RoadDistHBRWt, was similarly calculated within nested larger human population bins. The nal employed combination human population and road distance correction factor, PopRoadHBRWt, for a given milkweed observation was calculated as PopDenHBRWt × RoadDistHBRWt, and values were normalized such that PopRoadHBRWtN summed to the raw count per species (for further details, see Online Resource 1 Methods).
Corrected land cover data for each PopRoadHBRWtN weighted occurrence point was used to calculate the actual percentage of milkweed species per land cover, r i . We tabulated the 2016 NLCD land cover class raster cells falling within the MaxEnt core habitat to calculate the expected or available percent of land cover, p i . An omnibus chi-square goodness-of-t test was used to determine whether r signi cantly differed from p (P < 0.05), combining classes containing less than ve expected counts of milkweed occurrences. If the omnibus test was signi cant, we calculated 95% Bonferroni con dence intervals (BCIs) of r i for each land cover class and identi ed land cover classes where p i fell outside of the 95% BCI of r i (e.g., Neu et al. 1974;Byers et al. 1984;Manly et al. 2002;Kelly and Elle 2020).
To visualize differences in milkweed land cover electivity, we used the R electivity package (Quintans 2019) to calculate Vanderploeg and Scavia's (1979) relativized electivity index (E*; also referred to as Vanderploeg and Scavia's electivity index) for each land cover class: , where i = the i th land cover class, W i = (r i /p i )/(∑ r i /p i ) -1, r i is the proportion of the i th land cover class utilized, p i is the proportion of the i th land cover class available, and n is the number of land cover classes. The relativized electivity index incorporates several other electivity indices, including the forage ratio, r/p (Savage 1931;Ivlev 1961), and Chesson's (1978) alpha, W, which is also known as Vanderploeg and Scavia's (1979) selectivity coe cient and Manly's standardized selection index (Manly et al. 1972(Manly et al. , 2002Chesson 1978;Höner et al. 2010). Lechowicz (1982) recommended use of E* over several other electivity indices due to better approximation of a range from − 1 to 1, where − 1 indicates negative selectivity (avoidance), zero indicates random selectivity (neutral), and + 1 indicates positive selectivity (preference). For calculating combined species E i *, individual milkweed species E i * values were phenologically partitioned by multiplication with the species percentage occurrence for either the early or late 11-week period of monarch activity determined from the above milkweed phenology analysis. Phenologically partitioned values of E i * for the four species were then added for each bene t/risk category (for details on statistical analyses of E* in previous studies, the combination of E i * across species, construction of E i * category maps, and high E i * category kernel density estimates (KDEs), see Online Resource 1 Methods; for additional procedures, and scripts for R and ArcPython, see Online Resource 1 Table 3).

Niche modeling
The RSFSA results for niche model variable subset selection were used to determine that subsets of six to eight of the 95 variables were the optimal smallest sizes for producing higher performing MaxEnt models of the four milkweed species (Online Resource 1 Figs. 4-7, A-C). After generating MaxEnt models from thousands of random variable subsets of the selected subset sizes, the compared model selection criteria of AUC psa (pseudoabsence AUC used for most other niche modeling other than MaxEnt) and AICc bg (background AIC; Tracy et al. 2018) were both consistently effective in producing three replicates of 250 out of 3,000 models with signi cantly higher AUC psa and lower AICc bg than replicates of 250 random models out of 3,000 for all four milkweed species (Online

Core habitat distributions
The MaxEnt core habitat for A. a. ssp. capricornu primarily occupied the southeastern portion of the South-Central Semiarid Prairies from western Oklahoma to Central Texas (Fig. 2a). Asclepias viridis core habitat occurred along the eastern edge of the South-Central Semiarid Prairies from Kansas to Texas and included the southern Central Irregular Prairies, the eastern portion of the Ozark Ouachita Appalachian Forests (excluding the mountainous ecoregions), the western portion of the Southeastern USA Prairies, and the eastern Texas portion of the Texas Louisiana Coastal Plains (Fig. 2b). The MaxEnt core habitat of A. oenotheroides was concentrated along the southeastern edge of the South-Central Semiarid Prairies and northern portions of the Tamaulipas Texas Semiarid Plain. Asclepias oenotheroides core habitat included southern portions of the Texas Louisiana Coastal Plain, and portions of various Level II ecoregions in Mexico, including the Dry Gulf of Mexico Coastal Plains and Hills and Interior Depressions (Fig. 2c, Online Resource 1 Fig. 16). The core eastern habitat of A.
latifolia included the West-Central portion of the South-Central Semiarid Prairies from southwestern Kansas and southeastern Colorado to eastern New Mexico and northwest Texas, and portions of the northern Chihuahuan Desert (Fig. 2d).

Phenology relative to monarch activity
All four South-Central US milkweeds exhibited bimodal early and late season peaks in vegetative frequency. Asclepias a. ssp. capricornu and A. viridis had larger early-season peaks that corresponded and occurred more frequently during early season monarch activity from March to mid-May (Fig. 3). Asclepias oenotheroides had a relatively larger late season peak, and early and late season peaks were similar in magnitude for A. latifolia, and both these species occurred more frequently during late season monarch activity from mid-August to October (Fig. 4). The two seasonal peaks were associated with the rst and second early vegetative phenophases. Peak owering generally corresponded well with peaks in vegetative frequency, except for late season owering in A. latifolia, which peaked earlier than the overall vegetative peak (Figs. 3-4).
Land cover preference risk assessments The 2016 NLCD land cover classi cation had an overall accuracy of 43% at milkweed locations in the study area, and we used our manual land cover classi cations for calculating milkweed land cover utilization (r i ) (Online Resource 1 Table 17). The histogram bin ratio weighting sample bias adjustment effectively reduced the raw citizen science percentage utilization of roadside habitats for the four milkweed species by 24-55%, from 24-59% in the raw counts to 6-31% for the adjusted PopRoadHBRWtN counts (Fig. 5., Online Resource 1 Figs. 21-25). Roadside habitat utilization was signi cantly highest for A. latifolia, and A. asperula ssp. capricornu utilized roadsides to a signi cantly greater extent than A. viridis ( Fig. 5; Chi-squared tests for independence; P < 0.05).
The land cover classes with highest utilization by at least three of the four milkweed species were Grassland Herbaceous, followed by Developed, Open Space (Fig. 6a). These land covers generally dominated the most common general habitats utilized by the four milkweeds: Field/Meadow and Rural Road ROW.
Developed Medium and High intensity land covers also dominated within Rural Road ROW habitats (Online Resource 1 Figs. 21, 30). The land cover classes with the generally highest expected availability (p i ) in the milkweed core habitats included Grassland Herbaceous, Shrub/Scrub, and Cultivated Crops (Fig. 6b).
The two land cover classes of bene t (high E i *) for most milkweeds (except A. latifolia) were Developed-Open Space and Developed-Low Intensity (includes Developed, Medium and High intensity for A. oenotheroides; Fig. 6c). Cultivated Crops and Shrub/Scrub were prominent high-risk (low E i *) land cover for all four milkweeds (Fig. 6c) Intensity and Pasture Hay, and these were mapped as the most bene cial land covers in calculating the A. latifolia bene t KDE isopleths (Fig. 7d) and for calculating the milkweed species combination high E i *. Shrub/Scrub, Grassland Herbaceous, and Cultivated Crops land covers were all highly available and of high risk for A. latifolia (Fig. 6b-c), and they dominated most of its core habitat (Fig. 7d). Higher resolution views of the E i * land cover preference risk assessments revealed more detailed contrasts in land cover preference risk, such as for A. oenotheroides around urban areas of San Antonio and South Texas (e.g., Online Resource 1 Figs. 35-36).
Both the early and late season partitioned combination milkweed land cover preference risk assessment high E i * KDE surfaces and isopleths identi ed similar regional densities of milkweed high bene t land cover (Figs. 8-9). The combined high E i * 85% KDE isopleths (Figs. 8-9) resolved comparable major spring and fall interior and coastal corridors of land cover bene t for milkweed. The combined high E i * 85% KDE isopleth extended further into the Coastal Funnel in the fall compared to the spring (Fig. 9). The densest regions of milkweed land cover preference risk within the Central Funnel included core habitat for A. a. ssp. capricornu, A. viridis, and A. latifolia in western Oklahoma, A. latifolia habitat in West Texas, and A. oenotheroides habitat in South Texas (Figs. 7, 9). In the Coastal Funnel, the highest milkweed land cover preference risk was seen in South Texas A. oenotheroides habitat and central coastal Texas A. viridis and A. oenotheroides habitats (Figs. 7, 9). There was a potentially critical gap in milkweed habitat in South Texas within the Tamaulipas-Texas Semiarid Plain and southern Western Gulf Coastal Plain ecoregions, where much A. oenotheroides core habitat was at high land cover preference risk (Figs. 1, 7c, 8-9).

Niche modeling
The RSFSA revealed that subsets of six to eight of 95 variables were optimal for maximizing performance of MaxEnt models of the four milkweed species.
Tracy et al. (2019) and Kantola et al. (2019b) also found that six and eight variables was an optimal subset size for MaxEnt models using RSFSA with large sets of 80 and 119 variables, respectively. Consistent with previous niche modeling studies using RSFSA, we found signi cantly improved values of AUC psa and AICc bg for selected models over random models, but over tting was only sometimes improved (Tracy et al. 2018, 2019; Kantola et al. 2019b). We found that AUC psa criterion reduced over tting in more RSFSA replications than the AICc bg criterion, as found by Kantola et al. (2019b), but the AICc bg criterion can yield lower over tting with RSFSA in some cases (e.g., Tracy et al. 2018), or the two criteria may be equal in e cacy with regards to over tting (Tracy et al. 2019 Svancara et al. (2019). Our values of test AUC bgp were 5-8% lower than pseudoabsence AUC psa , and this was probably a result of lower false negative error rates that would arise in AUC psa since the pseudoabsence data were buffered 20-km from presence points. In contrast, AUC bgp considers as absences the background and presence points that can overlap the vicinity of true presence points to yield higher false negative error rates. MaxEnt default AUC bgp is more suitable for making comparisons between MaxEnt model studies, and AUC psa is best for comparisons of MaxEnt with other types of niche models, such as GLM.
Our feature selection results revealed the usefulness of incorporating additional climate variables in niche modeling beyond the typically employed 19 WorldClim indices (Bradie and Leung 2016), including the SuppClim and AET/PET indices, supporting our previous ndings (Tracy et al. 2018(Tracy et al. , 2019Kantola et al. 2019b). This is the rst study to nd that soil taxa are important variables for niche modeling of milkweeds, with soil taxa variables highly ranked in niche models of all four South-Central US species. Previous studies by Dilts et al. (2019) and Svancara et al. (2019) established that single edaphic variables were important in niche models of two studied western milkweed species, percent sand for A. a. ssp. asperula, and soil depth for A. speciosa, respectively. For example, we found that percent cover of usterts (Vertisols order) was the highest-ranking variable for niche modeling of A. oenotheroides, and the prevalence of usterts at a 1-km scale could be seen across the core habitat of A. oenotheroides in Texas (Online Resource 1, Table 4, Figs. 10a, 37).

Core habitat distributions and subspecies boundaries
Our models con rmed the importance of the South-Central US as core habitat for A. viridis, which included eastern Texas and Oklahoma and parts of western Arkansas and Louisiana but excluded most of Illinois and all of Indiana (Fig. 2b). These results contrast with the 0.5 probability contour for the MaxEnt distribution of A. viridis from Lemoine (2015) that extended past Illinois east to Indiana, but excluded most of Oklahoma, Texas, and Louisiana. The differences probably arise from our inclusion of a more representative number of A. viridis occurrence data to the West and South.
The 100% MaxEnt consensus model core habitat of A. a. ssp. capricornu generally resided well within the subspecies boundary, but it approached the boundary near the con uence of the Pecos River and Rio Grande in Texas (Fig. 2a). Our data supports the suggestion by Singhurst et al. (2015) that species status may be warranted for A. a. ssp. capricornu and A. a. ssp. asperula based upon clear differences in their morphology and distribution (see Woodson 1954). Further studies are needed to determine the distribution, frequency, and degree of hybridization in these subspeci c taxa, and our models indicated that it is important to include data from the southern Pecos River and Coahuila state.
Phenology relative to monarch activity Brower et al. (2018) found that early season monarch spring migratory activity in north central Florida coincided well with the vegetative phenology of the local milkweed monarch larval host, A. humistrata (pinewoods milkweed). Waterbury et al. (2019) found that phenology of Idaho and Washington A. speciosa was also well synchronized with early spring monarch arrival and oviposition. Similarly, we found that both early and late season monarch activity was generally well coordinated with early and late season peaks in milkweed vegetative phenology of the four studied milkweeds in the South-Central US. Asclepias latifolia represented an exception, with an early season vegetative peak occurring after monarch spring activity. However, the fall peak for A. latifolia corresponded well to late season monarch activity. Observed peaks in milkweed weekly vegetative phenology generally corresponded to seasonal peaks in owering phenology, except for the fall peak for A. latifolia, which had little owering.
The bimodal early and late season peaks in weekly vegetative frequency and the gaps between early and late season early vegetative phenophases for all four milkweeds (Figs. 3-4) may re ect late emergence of some individuals or resprouting following favorable late summer/early fall rainfall in certain years, which can be promoted by mid-summer disturbance. Calvert (1999) reported above average rainfall in the fall of 1996 produced abundant resprouting and peak October owering in all four studied species across Central Texas. We have also observed late September/early October owering of A. a. ssp. capricornu and A. oenotheroides in Central Texas following above average September 2018 rainfall without disturbance (Online Resource 1 Figs. 38-39). Nabhan et al. (2015) suggested that subspecies of A. asperula also may resprout following early summer damage. Late season resprouting of A. viridis has also been reported following summer rain (Dee and Palmer 2019) and summer disturbance, such as burning (Baum and Sharber 2012) or mowing (Baum and Mueller 2015;Dee and Baum 2019). Asclepias oenotheroides has the highest number of buds per root biomass of all four species, and A. latifolia is the only clonal species we studied (Pellissier et al. 2016). These more developed vegetative reproduction traits probably favor the observed relatively greater late season frequencies of A. oenotheroides and A. latifolia.
Land cover preference risk assessments The highest land cover preference risk was found with Cultivated Crops for all four milkweed species.  (Malcolm et al. 1993), was once abundant in cultivated elds (Yenish et al. 1997), where it was commonly utilized by monarch larvae, but where it now has been much reduced by increased herbicide use (Hartzler 2010;Pleasants and Oberhauser 2013;Zaya et al. 2017). Widespread herbicide use could also be lowering occurrence of A. oenotheroides and A. latifolia in Cultivated Crops in the South-Central US.
The Shrub/Scrub land cover class also represented high land cover preference risk for all four milkweeds, which could be related to lower soil moisture and anthropogenically lowered disturbance regimes in this land cover that fosters woody encroachment of open meadow milkweed habitats. Waterbury et al. (2019) found relatively low utilization of northwestern milkweeds for shrub/scrub, which they attributed to higher aridity. The Pasture/Hay land cover class represented high risk only for A. a. ssp. capricornu. The lower utilization of Pasture/Hay by A. a. ssp. capricornu may be related to a generally lower abundance of this species over the eastern portion of its range where Pasture/Hay is more abundant rather than to any lower suitability of Pasture/Hay for this milkweed compared to other species (Figs. 1, 2a, 6).
Grassland Herbaceous land cover was of high bene t for A. viridis and A. oenotheroides, no risk for A. a. ssp. capricornu, and high risk for A. latifolia (Fig. 6). Waterbury et al. (2019) found highest utilization of northwestern milkweeds for grassland/herbaceous cover. Asclepias latifolia favors grasslands, plains and prairies, trails, rocky canyon slopes, railways, and roadsides (Woodson 1954;Hart et al. 2000;Nabhan et al. 2015). Lower disturbance leading to less open areas in contemporary anthropogenically managed Grassland Herbaceous land covers may disfavor A. latifolia.
The Developed-Open Space and Developed-Low Intensity land covers, which commonly included roadsides and urban areas, were preferred by two to three of the four milkweed species, except for A. latifolia. Webb (2017) found that both A. viridis and A. a. ssp. capricornu occurred at higher frequency on roadsides compared to adjacent lands in Oklahoma. Asclepias viridis, along with A. syriaca and a variety of other milkweeds, are also common along roadsides in other states (e.g., Kaul et al. 1991;Kasten et al. 2016;Cariveau et al. 2020). In contrast, northwestern milkweeds had low utilization of developed land covers .
We endeavored to increase the accuracy of our land cover preference risk assessment by improving the quality of the input data sources through (1) using feature selected edapho-topo-climatic species distribution models to de ne core habitats, (2) using a manual land cover classi cation, (3) reducing sample bias in species occurrence surveys among the different land cover classes, and (4) minimizing temporal discrepancies between species occurrence data and the land cover classi cation. The overall accuracy of 43% for the 2016 NLCD classi cations at our milkweed occurrence sites was very low, justifying our decision to develop a manual land cover classi cation from Google Earth Pro and Google Map imagery. The low NLCD accuracy is probably due to the prevalence of milkweeds in edge habitats that are more prone to classi cation error, such as roadsides classi ed by NLCD as adjacent Grassland Herbaceous or Shrub/Scrub habitats rather than as Developed ROW land covers. The use of our novel histogram bin weighting sample bias reduction method lowered the frequency percentage of the four milkweed species in roadside ROW habitats by 24-55%, helping compensate for the natural bias of iNaturalist observations along roadsides (Online Resource 1 Figs. 22-25). However, there is still the potential for substantial error in the land cover risk assessment. For example, it is possible that our sample bias reduction method could be over-compensating for the potential bias of citizen science records near roads if a milkweed species has a greater a nity for disturbed habitats related to roadsides than the larger target group of nectar plants (see also Ranc et al. 2017).

Implications for monarch butter y conservation
The individual species and phenologically partitioned combined milkweed land cover preference risk assessments resolved major corridors of milkweed land cover bene t that may be especially important for providing monarch adult nectar and larval host plants for both the spring migration and the rst generation as well as the last summer generation and fall migration. The identi ed early and late season corridors of milkweed land cover bene t are similar and are mostly encompassed within the previously identi ed Central and Coastal monarch fall migratory funnels of Tracy et al. (2019). The spring monarch migration is dispersed, and spring migratory routes are more di cult to identify than fall migration routes (Solensky 2004), which have been clearly de ned (Calvert and Wagner 1999;Howard and Davis 2009;Tracy et al. 2019). Taylor (2019) described two major monarch spring migratory pathways through Mexico into Texas based upon spring rst sightings of Journey North, a primary interior pathway and minor coastal pathway. These spring migratory pathways correspond well with our fall migratory Central and Coastal funnels (Tracy et al. 2019), respectively (Fig. 9).
Large western portions of the Central Funnel consist of neutral to high-risk land cover for milkweed in both the early and late season, including areas in western Oklahoma and Central Texas where A. viridis and A. a. ssp. capricornu are more frequent during the early season and where A. latifolia is more frequent in West Texas during the late season ( Fig. 9). Increasing habitat for A. viridis and A. a. ssp. capricornu in western Oklahoma and west-central Texas would especially bene t the spring migration, while increasing A. latifolia habitat in West Texas should bene t late generation monarchs and fall migrants. A potential critical connectivity gap in milkweed land cover bene t was identi ed in A. oenotheroides core habitat . Providing additional A. oenotheroides milkweed habitat in South Texas may improve both spring and fall connectivity for the monarch migration.
We found from 6-31% of our milkweed occurrences were associated with roadsides in our study area, depending on the species. Management of roadside habitats is receiving increased focus for bene tting pollinator species in general (Hopwood et al. 2015;Phillips et al. 2020), and milkweeds and monarchs (Fischer et al. 2015;Kasten et al. 2016;Kaul and Wilsey 2019;Knight et al. 2019;Cariveau et al. 2019Cariveau et al. , 2020 Leone et al. (2019) reported that monarchs were more abundant at sites managed by prescribed re compared to grazing. Timing of mowing or burning in elds, roadsides, and urban open spaces can be especially important to maximize seasonal milkweed availability to monarchs for adult nectar resources and larval food plants (Baum and Sharber 2012;Baum and Mueller 2015;Fischer et al. 2015;Knight et al. 2019;Haan and Landis 2019a  Land cover class utilization/availability analysis (2014-2018) for Asclepias spp. in South-Central US within MaxEnt core habitats for (a) actual percent utilization of land cover classes (ri) by Asclepias spp. (n = locations weighted by human population and road distance, PopRoadHBRWtN); (b) expected percent availability (pi) for land cover classes, where values of pi with asterisks fall outside of the two-tailed 95% Bonferroni con dence intervals for ri (r signi cantly differed from p for all Asclepias spp.; Chi-squared Goodness of Fit Test, P < 0.05); and (c) land cover relativized electivity index (Ei*) for Asclepias spp., where land covers with signi cant differences noted in (b) are marked with colored arrows. Values of pi (b) with "x" had less than ve expected milkweed occurrences and were combined into the "Other" category, or Developed, High and Medium intensities were combined with Developed, Low Intensity when marked with "a", or Developed, High Intensity was combined with Developed, Medium Intensity when marked with "b" Asclepias spp. land cover preference risk assessments (2014-2018) according to land cover relativized electivity, Ei*, over MaxEnt core habitats (Fig. 4) within the South-Central US study area, including KDE isopleths for high milkweed land cover bene t Ei* and modi ed CEC (2005)  Combined Asclepias spp. land cover preference risk assessments (2014-2018) according to land cover relativized electivity index, Ei*, categories over MaxEnt core habitats and combined high Ei* 85% kernel density estimate (KDE) isopleths (Fig. 8)

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. TracyetalMilkweedsMSOnlineResource118Jul2021.docx