Distribution and phenology of monarch butterfly larvae and their milkweed hosts in the South Central US

South Central US milkweeds (Asclepias) 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). This study addresses multiple gaps in assessment of monarch conservation priorities for the South Central US through analyses of monarch larval host selectivity, phenology, and spatial density, as well as the phenology, niche modeled distribution, and land cover selectivity of important milkweed hosts. Results are synthesized to estimate seasonal milkweed resource areas. About 70% of monarch larval activity occurred from mid-March to mid-July (early season) and 30% from mid-August to late November (late season). Twenty-six wild milkweed (Apocynaceae) hosts were mapped, including four new records for North America. Important hosts included Asclepias a. ssp. capricornu, A. viridis, and A. oenotheroides, that were utilized more frequently during early season, and Asclepias latifolia, utilized more frequently during late season. Landscape host selectivity was positive for A. viridis and A. a. ssp. capricornu in late and early seasons, respectively, and negative for A. oenotheroides in late season. Milkweed land cover selectivity was positive for Developed-Open Space and Grassland Herbaceous, and negative for Cultivated Crops and Shrub/Scrub. Seasonal milkweed resource areas and larval spatial densities resolved interior and coastal corridors providing functional connectivity for monarch spring and fall migrations. A potential gap in milkweed resource areas was identified in South Texas. The novel merging of milkweed niche models with larval phenology, host selectivity, milkweed phenology, and land cover selectivity informs conservation assessment.


Introduction
The monarch butterfly (Danaus plexippus) recently became a candidate for US listing as federally endangered or threatened (USDI FWS 2020). Eastern migratory monarchs are unique among insects in flying up to 3000 km in their fall migration from breeding sites as far north as southern Canada to overwintering sites in Central Mexico (Urquhart and Urquhart 1978), with the same butterflies returning to the US in the spring. Spring remigrants lay eggs in the South Central US from Texas to Oklahoma on milkweeds (Asclepias spp.), their primary larval host plants (Malcolm et al. 1993). Batalden et al. (2007) and Castañeda et al. (2019) identified the South Central US as the primary region of activity for first generation monarchs from March to April. Therefore, milkweed habitat in the South Central US (Texas, Oklahoma, Louisiana, and Arkansas) is critical for development of the first generation in the spring, which forms the basis for subsequent generations that migrate and reproduce further north (Flockhart et al. 2013;Oberhauser et al. 2017). Milkweeds are also preferred nectar plants of adult monarchs (USDA NRCS 2015a, b, 2018 and can be important for fueling the spring and fall migrations (Brower et al. 2006;Hobson et al. 2020). Because of the importance of South Central US milkweed habitat for maintaining monarch functional connectivity along the migratory corridors, spatially explicit niche models (Schwartz 2012) and habitat selectivity assessments (Johnson 1980) in this region are needed to inform monarch conservation planning efforts (e.g., Moreno-Sanchez et al. 2019). The significant long-term decline in the eastern monarch population that overwinters in Central Mexico (Semmens et al. 2016) has been primarily attributed to the loss of milkweed habitat in the Midwest US for second and third generation larvae Pleasants 2017;Thogmartin et al. 2017;Stenoien et al. 2018;Taylor et al. 2020), adult mortality within the Central Mexico overwintering grounds (Brower et al. 2004(Brower et al. , 2012, breeding season climate change (Zylstra et al. 2021), or more of a combination of these factors (Wilcox et al. 2019). Demographic models of Oberhauser et al. (2017) indicate that increasing milkweed and nectar availability for monarchs in both the northern and southern parts of the breeding range may be the most efficient strategy for promoting monarch population growth. Few studies have focused on monarch conservation with respect to the South Central US (but see Tracy et al. 2019;Kantola et al. 2019a). We investigate South Central US monarch larval and milkweed distribution and phenology, advancing several new applications of monarch community science data (formerly "citizen science data") to those reviewed by Ries and Oberhauser (2015).
Analysis of monarch host selectivity at a landscape level can be valuable in understanding the relative importance of host milkweeds in the South Central US. Few studies have analyzed insect population host selectivity at a landscape level. Lechowicz and Jobin (1983) and Mauffette et al. (1983) investigated landscape host selectivity of gypsy moth larvae for tree species in Quebec forests using the relativized electivity index (E*, Vanderploeg and Scavia 1979). Martinson et al. (2016) also utilized E* in evaluating landscape host selectivity of brown marmorated stink bugs for woody plant species in nurseries. We expand the use of E* and estimate monarch larval landscape host selectivity for milkweeds in the South Central US, uniquely employing community science data to reveal important host associations. We also novelly employ community science larval occurrence data to assess monarch larval phenology and spatial density distribution across the South Central US, incorporating these data into seasonal landscape host selectivity indices. Understanding the important milkweed host species and peak 1 3 seasonal periods of monarch larval activity is fundamental for guiding an evaluation of the spatial distribution and phenology of important milkweed resources for monarch butterflies in the South Central US.
An understanding of the influence of land use/land cover (LULC) classes on milkweed habitat is critical for estimating the distribution of milkweed resources for monarchs in the South Central US. Asclepias syriaca favors agricultural landscapes in Canada (Pitman et al. 2018) and roadsides in Iowa (Kaul and Wilsey 2019), but it has significantly declined in the Midwest over the last century in conjunction with agricultural intensification (Hartzler 2010;Boyle et al. 2019) related to increased use of glyphosate herbicides (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. We uniquely employ E* to assess milkweed land cover selectivity in the South Central US with community science data. To our knowledge, Ogawa and Mototani (1991) performed the only other study that applied utilization/availability analysis to assess plant land cover selectivity (see Supplement Introduction on utilization/availability analysis).
Niche models for North America milkweed species have been essential to investigate the actual and potential geographic distributions of adult and larval monarch food resources under current climate in the western US (Dilts et al. 2019), and both current and future climate scenarios for the northeastern US (Lemoine 2015) and northwestern US . Kass et al. (2020) incorporated milkweed niche models into modeling the monarch fall migration from South Texas to Mexico. Except for range-wide MaxEnt models developed for A. viridis (Lemoine 2015), range-wide models for South Central US milkweed species are lacking, and we develop them for this study. We take a hierarchical approach (see Wilting et al. 2010 and Supplement Introduction) to combine milkweed niche model output with results from phenological and landscape resource selectivity analyses to estimate seasonal milkweed resource areas across the South Central US.
The objectives of this study were to (1) analyze monarch larval landscape host selectivity and larval phenology to identify important focal milkweed hosts and the major seasonal periods of monarch larval activity in the South Central US; (2) delineate seasonal monarch larval spatial density distributions; (3) develop MaxEnt niche models for defining core habitats of focal milkweed species; (4) compare seasonal phenology of focal milkweeds; (5) assess milkweed land cover selectivity; and (6) estimate seasonal milkweed resource areas for monarch butterflies through combining the results from all analyses. Seasonal milkweed resource areas are compared with monarch larval spatial densities, and results are discussed in the context of monarch conservation.

Study area
The target study area comprised the South Central US region, including Arkansas and the majority of Oklahoma, Texas, and Louisiana, which represents the principal area for the development of first and last generation monarch larvae ( Fig. 1; see also Fig. S1). The western edge of the study area consisted of the western distribution limit for A. a. ssp. capricornu in West Texas and Oklahoma (Woodson 1954) (see Occurrence data, below). The eastern boundary consisted of the eastern borders of Arkansas and northern Louisiana, and the Mississippi River in southern Louisiana (Fig. 1). The background evaluation extents for niche modeling each of the focal 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.

Occurrence data
A variety of data sources were used in six primary analyses for developing seasonal milkweed resource area maps for the South Central US (Fig. 2, items 1-6). Monarch larval instar associations with wild (unplanted) native hosts by locality and date for the South Central US study were compiled for 2011-2020 from iNaturalist (2021a) and for 2011-May 2017 from Journey North (2017) and Monarch Larva Monitoring Project (2017). Larval data from extensive seasonal focused surveys (e.g., Berman and Baum, unpublished data) were not included to reduce seasonal survey bias. To reduce spatiotemporal bias, the larval data were spatially filtered by host species to 10 km (Boria et al. 2014) leaving 341 occurrences for the landscape host selectivity analysis (Fig. 2, item 1). Spatio-temporal filtering by week/year within 30 km per host left 547 occurrences for the larval phenology analysis by hosts (Fig. 2, item 2). Both sets of larval data were then weighted to reduce sample bias in higher human population density areas (see Phenology of monarch larvae, below). Results from the first two analyses were used in selecting the focal milkweed hosts Fig. 1 South Central US study area with Danaus plexippus larvae and pupae (2011-2020) occurrences on wild native Apocynaceae (n = 341 spatio-temporally thinned locations used in analyses), with selected additional occurrences (n = 19, data outside of analyzed spatio-temporal boundaries, see Table S1), and D. plexippus percent landscape utilization of host species (inset chart), with fall D. plexippus migratory funnels from Tracy et al. (2019) (see Table S4 for host and distribution records) 1 3

Fig. 2
Flowchart of various analysis output items (circled digits, e.g., ❶) for estimating the distribution of seasonal milkweed resource areas 1 3 and defining periods of monarch larval seasonal activity for the four remaining analyses (Fig. 2, items 3-6). The 10 km thinned larval distribution data (Fig. 2, item 1) were also used to model seasonal larval spatial density distributions for early (weeks 11-28, 232 occurrences) and late (weeks 32-47; 195 occurrences) seasons (Fig. 2, item 3). For assessment of novel monarch larval native host distribution records in North America, we surveyed the literature to update host distribution records following Malcolm and Brower (1986), Robinson et al. (2002), and Beccaloni et al. (2008). For developing niche models and defining core habitat for selected focal milkweed hosts (Fig. 2, item 4), historical (1834-2018) occurrence data throughout the entire North American ranges for the focal milkweed species were obtained from a variety of sources. Milkweed occurrence sources included iNaturalist (2021a), SEINet (2017) ]. An effort was made to obtain occurrence data representative of the known county distributions of milkweed species (Kartesz 2015), such as through georeferencing SEINet, vPlants, and OVPD records in various counties from which we lacked specific locations. From 388 to 979 occurrence points per species were available for modeling after spatial thinning to 10 km to reduce spatial autocorrelation (Boria et al. 2014). We partitioned data for A. a. ssp. capricornu to the East from A. a. ssp. asperula to the West to improve niche model quality (e.g., Gonzalez et al. 2011; for details on the subspecies boundaries, see Supplementary Information). To reduce spatiotemporal bias for focal milkweed species phenology analyses (Fig. 2, item 5), 2011-2020 iNaturalist milkweed records over their identified core habitats in the South Central US were spatiotemporally filtered by week/year within 30 km per host to leave 272 to 520 occurrences. In the milkweed land cover selectivity analyses (Fig. 2, item 6), land covers were manually determined for 2014-2018 iNaturalist focal milkweed occurrence points (within 2 years of 2016 land cover data) with at least ca. 30-m resolution within MaxEnt modeled core habitat for the South Central US study. Milkweed land cover class data were spatially thinned to 1 km to leave from 67 to 586 occurrence points, which were sample bias adjusted for human population density and road distance (see below Phenology of monarch larvae). Milkweed land cover classes were assigned based on the 2016 National Land Cover Database (NLCD) legend (Jin et al. 2019;MLRC 2021a, b) within a 15 m radius of the location and within the context of the adjacent surroundings according to historical imagery for 2014-2018 from Google Earth Pro (2021) and Google Maps (2021) Street View. In addition, the general habitat of the milkweed location was manually classified using the imagery into road rights-of-way (ROW) versus non-road habitats (for further details on manual land cover classifications using imagery and further general habitat subdivisions, see Supplementary Information; for all occurrence data Excel files, see Table S1).

Phenology of monarch larvae and milkweeds
The phenology of both monarch larvae and milkweed host plants was assessed for the 2011-2020 period in the South Central US (Fig. 2, items 2 and 5). Spatiotemporally thinned occurrence data were used to develop weekly stacked phenophase histograms and phenoperiod Gantt charts. Monarch larval phenophase histograms were visually inspected 1 3 to demarcate between early and late season activity periods. The early and late season larval and milkweed frequency data were fitted to curves with highest R 2 using CurveExpert software (Hyams 2010) (for further details on determination of monarch larval hosts and phenophases of monarch larvae and milkweeds from community science data images and statistical analysis, see Supplementary Information; for excel histogram charts and R script for analyses, see Tables S1 and S3, respectively).
The estimated relative proportions of community science observations of monarch larvae on various milkweeds predominating near and away from human population centers are sensitive to errors from sample bias (e.g., Geldmann et al. 2016;Zhang and Zhu 2018). Larval-host association occurrences were weighted to reduce sample bias using a novel target group-based Histogram Bin Ratio Weighting (HBRWt) method (sample bias method hereafter). 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 community 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., human population density) 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 for correcting sample bias that leads to more conservative estimates of milkweed host associations within human population dense areas. The HBRWt sample bias method involved first constructing histograms of human population density for two sets of data in the study area: (1) 10,000 unthinned random points; and (2) up to 10,000 unthinned target group points, where the target group was 2011-2020 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). Values of PopDenHBRWt were normalized to yield PopDenHBRWtN which summed to the raw count per species (for further details, see Supplementary Information).

Landscape resource selectivity analyses
To assess monarch landscape host selectivity and milkweed land cover selectivity (Fig. 2, items 1 and 6), we employed a utilization/availability analysis that incorporated Vanderploeg and Scavia's (1979) relativized electivity index, E*. The R electivity package (Quintans 2019) was used to calculate E* for each larval host or land cover class: is the proportion of the ith host or land cover class utilized, p i is the proportion of the ith host or land cover class available, and n is number of hosts or land cover classes. 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, zero indicates random selectivity, and + 1 indicates positive selectivity. The PopDenHBRWtN weighted larval occurrence points were used to estimate the actual percentage of larval occurrences per milkweed host, r i . In addition to using PopDenHBRWt, we also employed a road distance HBRWt correction factor, RoadDistHBRWt, for adjusting milkweed land cover occurrences (only available for the higher resolution land cover data) that was calculated similarly to PopDenHBRWt within nested larger human population bins. A combined human population and road distance correction factor, PopRoadHBRWt, for a given milkweed land cover observation was calculated as PopDenHBRWt × RoadDistHBRWt, and values were normalized to sum to the raw count per species (for further details, see Supplementary Information, Landscape resource selectivity analyses). The manually determined land cover data for each PopRoadHBRWtN weighted milkweed land cover occurrence point were used to estimate the actual percentage of milkweed species occurrences per land cover, ri. Host plant occurrence records (iNaturalist 2021a) for 2011-2020 were combined with larval host record data, thinned by 10 km per species, and weighted by PopDenHBRWt to estimate the percent available hosts, p i . For the expected or available percent of land cover, p i , we tabulated the 2016 NLCD land cover class raster cells falling within the MaxEnt core habitat (see Results for p i sample sizes). The manually determined occurrence point land cover classifications for milkweed occurrences used for r i were also compared with classifications derived from the NLCD 2016 spatial database. For the monarch landscape host selectivity analysis, data for r i and p i were partitioned into three periods prior to thinning for seasonal analyses: year-round, early season, and late season (see Phenology above). An omnibus chi-square goodness-of-fit test was used to determine whether r i significantly differed from p i (P < 0.05), combining classes containing less than five expected counts of host or milkweed land cover occurrences. If the omnibus test was significant, we calculated 95% Bonferroni confidence intervals (BCIs) of r i for each host or land cover class, and we reported host and 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). Habitat utilization (e.g., roadside versus non-road habitats) was also analyzed for milkweed land cover observations using Chi-squared tests for independence.

Monarch larval spatial density distribution
Kernel density estimation of species occurrences can be useful in representing species spatial distributions (Pennay et al. 2011), especially for migratory species such as the monarch butterfly (Tracy et al. 2019). The early and late season spatial densities of monarch larval observations were modeled through developing 1 km resolution kernel density estimation (KDE) surfaces (geodesic method, expected counts) in ArcGIS using the 10 km thinned and sample bias adjusted larval occurrence data (see above) ( Fig. 2, item 3). The monarch larval KDE 85% isopleth was used to visualize most of the monarch larval density distribution.

Milkweed niche modeling
An initial set of 95 environmental variables at 1 km resolution was used in developing MaxEnt niche models (Fig. 2, item 4), 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) (Table S2). Anthropogenic and land cover variables were intentionally excluded from the niche models to facilitate a more independent assessment of species affiliation with land cover through the milkweed land cover selectivity analysis (see above) (for justification of initial variable selection, see Supplementary Information). MaxEnt (Phillips et al. 2006) version 3.3.3 models were developed from milkweed locations using the R-software (R Core Team 2017) dismo package (Hijmans et al. 2011). We generally followed methods described in Tracy et al. (2018) for generation of 10,000 background and 10,000 pseudoabsence points and calculation of the performance metrics of Area under the Curve (AUC), True Skill Statistic (TSS), and overfitting (for details, see Supplementary Information).
We utilized the random subset feature selection algorithm (RSFSA) to evaluate performance of thousands of MaxEnt milkweed species 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(Tracy et al. , 2019Kantola et al. 2019b) (for further details, see Supplement, Methods and Tracy et al. 2018Tracy et al. , 2019. A total of 48 feature selected binary calibrated MaxEnt models were projected for mapping, 12 for each of the four species. These 12 selected models per species were assembled into feature subset ensembles, with the area of 100% model frequency consensus used to represent milkweed core habitats (e.g., Porfirio et al. 2014) (for further details on niche model calibration, consensus, and variable ranking, see Supplementary Information and Tables S12-S15).

Seasonal milkweed resource areas
Results from the above analyses (Fig. 2, items 1-6) were synthesized in the form of various seasonal indices representing milkweed species habitat quality that were applied to milkweed core habitat areas for raster visualization of seasonal milkweed resource areas (Fig. 2, items 7-10). Output from the monarch larval host selectivity and phenology analyses and the milkweed core habitat phenology analysis (Fig. 2, items 1-2, 4-5) were first summarized in the form of a seasonal landscape host selectivity index (SLHSeI) for each milkweed host (Fig. 2, item 7), calculated as SLHSeI = E i * SNH_T × %STF x %SHF, where E i * SH_T (Fig. 2, item 1) is the seasonal E i * value representing selectivity for a host by the target organism that has been integer normalized from zero to 10,000 to facilitate multiplication and raster assignments, %STF (Fig. 2, item 2) is percent seasonal target organism (monarch larva) frequency, and %SHF (Fig. 2, items 5) is the seasonal host (milkweed) frequency in the core habitat (Fig. 2, item 4). The SLHSeI values were then combined with the milkweed land cover selectivity results, to produce the land cover/seasonal landscape host selectivity index (LCSLHSeI) for each focal milkweed species (Fig. 2, item 8), calculated as LCSLHSeI = SLHSeI × E i * LC_H , where SLHSeI is defined above and E i * LC_H is the E* i value representing land cover selection by the host that has been integer normalized from zero to 10,000 to facilitate multiplication and raster value assignments (index calculation tables are referenced in the Results). Values of LCSLHSeI were assigned to each land cover class pixel within core habitats for each focal milkweed species using the reclassify tool in ArcGIS to create three separate rasters per milkweed species per season (three early season and three late season rasters), one raster for each land cover class selectivity category of positive, random, or negative (Fig. 2, item 8). The resulting three categorized milkweed LCSLHSeI rasters per season were added by category across all species using Arc GIS cell statistics summation (with "no data" ignored) to create combined seasonal LCSLHSeI rasters (Fig. 2, item 9). Kernel density estimate surfaces were also developed in Arc GIS from early and late season combined positive selectivity category LCSLHSeI (milkweed resource indices), including 85% KDE isopleths to represent seasonal milkweed resource areas (Fig. 2, item 10). Results were displayed with 85% KDE isopleths of seasonal larval spatial density (Fig. 2, item 3) (see Supplementary Information; for details on kernel density estimate analysis and additional summary procedures, and scripts for R and ArcPython, see Table S3).

Monarch landscape host selectivity
Twenty-six wild native monarch larval hosts were identified across the South Central US (non-native ornamental A. curassavica was not included in this study). These included four new host records, A. linearis, A. arenaria, A. brachystephana, and Cynanchum unifarium (Fig. 1), according to our updated North American monarch larval host distribution records (Table S4). In addition, 18 new state host distribution records were found among 13 previously known host species, including A. asperula ssp. capricornu, A. oenotheroides, A. tuberosa, A. syriaca, A. variegata, A. purpurascens, A. subverticillata, A. incarnata, A. speciosa, A. sullivantii, A. engelmanniana, A. quadrifolia, and C. laeve (Fig. 1, Table S4). The top six hosts utilized annually (in order of utilization) were A. viridis, A. asperula ssp. capricornu, A. oenotheroides, C. laeve, A. latifolia, and A. tuberosa. The four top Asclepias spp. all had higher availability than C. laeve, and they were chosen as focal species in this study (Figs. 1, 3A, B). The top five available hosts across the region annually were A. viridis, A. asperula ssp. capricornu, A. oenotheroides, A. tuberosa, and A. viridiflora (Fig. 3B). We observed seven other milkweed larval host records in South Central US states from species from which monarch larvae have previously been reported, including A. perennis, A. texana, A. amplexicaulis, A. hirtella, A. longifolia, A. viridiflora, and Matelea reticulata (Fig. 1, Table S4). Only three previously recorded monarch larval hosts for the South Central US were not seen in data from this study, Funastrum cynanchoides, F. crispum, and A. obovata (Table S4). Annual positive monarch landscape host selectivity was identified for A. viridis and A. asperula ssp. capricornu, while negative selectivity was identified for A. tuberosa, A. viridiflora, A. amplexicaulis, A. engelmanniana, and Matelea reticulata (Fig. 3C). Several differences were noted in seasonal landscape host selectivity, such as negative selectivity for A. oenotheroides in the late season, and negative selectivity for A. quadrifolia and A. variegata in the early season (Fig. 3C) (see below for seasons).

Phenology of monarch larvae and milkweeds
Early and late season peaks of monarch larval activity in the South Central US were centered around late April and early October, respectively. Early and late season periods of monarch larval activity extended from mid-March to mid-July (weeks 11-28) and mid-August to late November (weeks 32-47), respectively, with a general between season cutoff around early August (weeks 31-32) (Fig. 4). Larval frequency was significantly higher in the early season for all the top four utilized Asclepias hosts except A. latifolia, for which larval activity was significantly highest in the late season (Fig. 4B).
All four focal South Central US milkweeds exhibited bimodal early and late season peaks in vegetative frequency. Asclepias viridis and A. a. ssp. capricornu had larger early season peaks that corresponded well with early season monarch activity from mid-March to mid-July (Figs. 4, 5A, B). The relative frequency of A. oenotheroides and A. latifolia did not significantly differ between early and late seasons (Fig. 5C, D). The two seasonal peaks in focal milkweeds were associated with the first and second early vegetative phenophases. Peak flowering generally corresponded well with peaks in vegetative frequency ( Fig. 5A-C), except for late season flowering in A. latifolia, which peaked earlier than the overall vegetative peak (Fig. 5D).  PopDemHBRWtN); B expected percent availability (p i ) for milkweed hosts, where values of p i with asterisks fall outside of the two-tailed 95% Bonferroni confidence intervals for r i (r significantly differed from p for all seasons; Chi-squared Goodness of Fit Test, P < 0.05); and C larval milkweed host relativized electivity index (E i *), where milkweed hosts with significant differences noted in (B) are marked with colored arrows. Values of p i (B) with "x" had less than five expected milkweed occurrences and were combined into the "Other" category, or various host "Grp" categories ( a A. quadrifolia species present as separate category and not in group for early season; see Table S5 for groupings per period)

Monarch larval spatial density distribution
Both early and late season community science monarch larval host records from 2011 to 2020 were heavily biased towards population centers, such as San Antonio/Austin, Dallas/ Fort Worth, Houston, and Oklahoma City (Fig. 6). The highest larval densities occurred along the eastern edge of the Central Funnel and central portion of the Coastal Funnel. Primary ecoregions with higher larval densities were the eastern edge of the South Central Semiarid Prairies, the central portions of the East Central Texas Plains and Texas Louisiana Coastal Plain, and a western portion of the Ozark Ouachita Appalachian Forests-West (Figs. 6, S1). Late season monarch larval densities were much lower than those in the early season, and the densities (85% KDE isopleths) were more widely distributed to the West, North, and East (Fig. 6).

Milkweed niche modeling
Milkweed niche model performance was high for all four focal milkweeds in terms of test pseudoabsence AUC (AUC psa ; AUC used in feature selection) (0.86-0.97), and moderate to high in terms of background AUC (AUC bgp ; default MaxEnt AUC) (0.81-0.9) (Table S11). Model AUC psa overfitting (AUC psa_diff ) was similar across the four species, ranging from 0.002 to 0.011, giving an average value of 0.008, which is a little lower than the average value for AUC bgp overfitting (AUC bgp_diff ) of 0.013 (Table S11). The random subset feature selection algorithm analysis revealed 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 (Figs. S11-S14, A-C) (see Supplement Results for details on variable importance and feature selection).
Asclepias viridis MaxEnt core habitat occurred along the eastern edge of the South Central Semiarid Prairies from Kansas to Texas and included the southern Central Irregular Plains, the eastern portion of the Ozark Ouachita Appalachian Forests (excluding the mountainous ecoregions), the western portion of the Southeastern USA Plains, and the eastern Texas portion of the Texas Louisiana Coastal Plain (Fig. 7A). 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. 7B). Asclepias oenotheroides core habitat was concentrated along the southeastern edge of the South Central Semiarid Prairies and northern portions of the Tamaulipas Texas Semiarid Plain, extending southwards into Mexico along the East coast (Fig. 7C, S23). 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. 7D).

Fig. 4
Danaus plexippus larvae and pupae annual (2011-2020) weekly stacked phenophase histograms and phenoperiod Gantt charts (= legend) for South Central US study area (n = spatio-temporally thinned locations), with percent of total area under fitted nonlinear curves for the monarch immature activity periods of mid-March to mid-July and mid-August to late November for: A D. plexippus larval and pupal stages, and B native Apocynaceae host plants (for larger images, fitted curve formulae, and R 2 values, see Figs. S5, S6). The two-tailed 95% Bonferroni confidence intervals for early and late season frequency percentages with asterisks fall outside of the expected random frequency of 50%, assuming they were equal (omnibus Chisquared Goodness of Fit Tests, P < 0.05)

Milkweed land cover selectivity
The 2016 NLCD land cover classification had an overall accuracy of 43% at milkweed locations in the study area, and we used our manual occurrence point land cover classifications for calculating milkweed land cover utilization (r i ) ( Table S24). The sample bias adjustment effectively reduced the community science percentage proportions of milkweeds in roadside habitats downward from 24-59% in the raw counts to 6-31% for the adjusted PopRoadHBRWtN counts (Figs. S24-S29). Roadside habitat utilization was significantly higher for A. latifolia, and A. asperula ssp. capricornu utilized roadsides to a significantly greater extent than A. viridis ( Fig. S24; Chi-squared tests for independence; P < 0.05).
The land cover classes with highest utilization by at least three of the four focal milkweed species were Grassland Herbaceous, followed by Developed Low, Open Space (Fig. 8A). Field/Meadow and Rural Road ROW land covers generally dominated the most common habitats utilized by the four milkweeds (Figs. S30D, S32D, S32D, S33D). Developed Medium and High Intensity land covers also dominated within Rural Road ROW habitats (Fig. S34). The land cover classes with the generally highest availability (p i ) in the milkweed core habitats included Grassland Herbaceous, Shrub/Scrub, and Cultivated  (2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) weekly stacked phenophase histograms and phenoperiod Gantt charts (= legend) for MaxEnt core habitats within South Central US study area (n = spatio-temporally thinned locations), with percent of total area under fitted nonlinear curves for the Danaus plexippus activity periods of mid-March to mid-July and mid-August to late November for: A A. viridis, B A. asperula ssp. capricornu, C A. oenotheroides, and D A. latifolia (for larger images, fitted curve formulae, and R 2 values, see Figs. S7-S10). The two-tailed 95% Bonferroni confidence intervals for early and late season frequency percentages with asterisks fall outside of the expected random frequency of 50%, assuming they were equal (omnibus Chi-squared Goodness of Fit Tests, P < 0.05) Crops (Fig. 8B). The two land cover classes of positive selectivity for most focal milkweeds were Developed-Open Space and Developed-Low Intensity. Cultivated Crops and Shrub/Scrub were prominent negative selectivity land cover for all four focal milkweeds (Fig. 8C).

Seasonal milkweed resource areas
Both the early and late season milkweed resource areas occurred over similar regions along the eastern half of the Central Funnel and a central portion of the Coastal Funnel (see Tables S6-S10 for values used in calculation of milkweed resource indices). The milkweed resource area densities were much higher in the spring (Fig. 9). For example, the late season milkweed resource area 5% KDE isopleth, which matches the density value of the early season milkweed resource area 85% KDE isopleth, only occupied two small areas in the east central portion of the Texas Central Funnel (Fig. 9B). The late season milkweed resource areas extended a little further west and east compared to the spring (Fig. 9). The early season larval 85% KDE isopleth regions (larval density areas, hereafter) mostly fell within the milkweed resource areas in the eastern Central Funnel and a portion of the central Coastal Funnel. The late season larval density areas overlapped similar late season milkweed resource areas, but portions extended further westwards and eastwards outside of the main milkweed resource area regions. In both early and late season milkweed resource areas, there was a potentially critical gap in South Texas, where there were much negative selectivity milkweed resource index areas in A. oenotheroides core habitat ( Fig. 7C; southernmost red areas in Fig. 9).  PopRoadHBRWtN); B expected percent availability (p i ) for land cover classes, where values of p i with asterisks fall outside of the two-tailed 95% Bonferroni confidence intervals for r i (r significantly differed from p for all Asclepias spp.; Chi-squared Goodness of Fit Test, P < 0.05); and C land cover relativized electivity index (E i *) for Asclepias spp., where land covers with significant differences noted in (B) are marked with colored arrows. Marked values of p i (B) were combined as follows: "x" had less than five expected milkweed occurrences and were combined into the "Other" category, "a" had Developed, High and Medium Intensities combined with Developed, Low Intensity, "b" had Developed, High Intensity combined with Developed, Medium Intensity 1 3

Monarch landscape host selectivity
Community science data on monarch larval wild host associations proved valuable in revealing regionally important and novel host associations. The four new native host records bring the total to 58 host species for North America, of which about half (29) are native in the South Central States, including the 26 species found in this study (Table S4). We selected A. viridis, A. a. ssp. capricornu, A. oenotheroides, and A. latifolia as the most important focal milkweeds for this study as they had high landscape host utilization and availability. These four species were the primary late season monarch larval hosts reported in Central Texas by Calvert (1999). Lynch and Martin (1993)  Only two milkweeds had significant positive landscape host selectivity, A. viridis in the late season and A. a. ssp. capricornu in the early season. Significant negative host selectivity for A. oenotheroides in the late season indicates it might be a less preferred or suitable host than A. viridis during that time. A variety of other ecological interactions could also account for lower selectivity for A. oenotheroides, including higher predation or parasitism (e.g., Kursar et al. 2006) that could be associated with habitats favored by A. oenotheroides, including habitats where resprouting following disturbance tends to occur. Predation of monarch larvae can be higher in grassland (Myers et al. 2019) and urban habitats (Baker and Potter 2020), and A. oenotheroides exhibited significant positive selectivity for the Developed, Low Intensity land cover class that was not seen in A. viridis (Fig. 8C). Asclepias tuberosa had negative selectivity for year-round and seasonal periods, which may be partly related to lower monarch field ovipositional preference for A. tuberosa compared to several other Asclepias spp. (Pocius et al. 2018). Although A. viridiflora was among the top five available hosts for the region, it had negative year-round and late season host selectivity, and further studies are needed on monarch ovipositional preference and larval suitability for this species. Cynanchum laeve was among the top five utilized hosts in the South Central US, but it had lower availability compared to the four focal Asclepias spp. Monarch larval survival on C. laeve is similar to that on other highly suitable Asclepias hosts (Yeargen and Allard 2005; Pocius et al. 2017). Cynanchum laeve appears to not only be a more important monarch host plant than generally realized for the East Central US (Yeargen and Allard 2005), but also for the South Central US.
There is probably still significant bias towards monarch larval hosts encountered in proximity to urban areas in our sample bias adjusted community science data. Regional community science groups such as the "Texas Milkweeds and Monarchs" iNaturalist project (iNaturalist 2021b) collected valuable data for our study but can also be a potential source of spatial bias. Consequently, results from our landscape host selectivity analysis should be interpreted with caution and compared with data from field studies with random sampling protocols. Monarch landscape host selectivity analyses can reveal important host associations in other areas with monarch community science data, including Mexico.

Phenology of monarch larvae and milkweeds
The community science phenology data are apparently biased towards more visible larger fifth instar monarch larvae (Fig. 4), and possibly towards larger, and more showy flowering milkweed host stages or species (Fig. 5). These observer biases may have skewed seasonal peak frequencies in favor of weeks with higher counts of more apparent stages. Consequently, caution is warranted in interpreting the phenology results, but they can indicate primary seasonal activity periods. Lynch and Martin (1993) found that April-May first generation monarch larval activity in Texas and Louisiana was best synchronized with the presence of A. viridis and A. asperula ssp. capricornu. They suggested that a second, and possibly third, generation of monarch larvae may be produced past May in the northern South Central US. Second generation larvae may dominate our June to mid-July early season monarch larval records. Early season monarch spring migratory activity was also well synchronized with local Asclepias spp. phenology in North Central Florida (Brower et al. 2018) and Idaho and Washington . Similarly, we found that both early and late season peaks in monarch larval activity of late May and early October generally agreed well with early and late season peaks in milkweed vegetative phenology of the four focal milkweeds in the South Central US. Asclepias latifolia represented an exception, with an early season vegetative peak occurring in early July, much later than peak late April monarch spring activity. However, the late September peak for A. latifolia frequency approximately corresponded to late season early October peak in monarch larval activity. Late season monarch larval activity in the South Central US can represent around 30% of the total seasonal monarch larval frequency (Fig. 4). Further investigation is needed to determine the degree to which late season larvae represent fourth or later generations from midsummer migrants originating north of the South Central US or from resident populations in South Central US urban areas , which potentially breed on ornamental A. curassavica. The phenology of monarch larvae on ornamental A. curassavica warrants investigation. Finer spatio-temporal partitioning of the monarch larval and milkweed phenology data can be used to investigate regional variation in phenology and changes in phenology over time.
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 (Fig. 5) may reflect 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 flowering in all four studied species across Central Texas. We have also observed late September/early October flowering of A. a. ssp. capricornu and A. oenotheroides in Central Texas following above average September 2018 rainfall without disturbance (Figs. S36, S37). 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.

Monarch larval spatial density distribution
Patterns of monarch larval density were generally similar between the early and late season and mainly centered around host records on A. viridis and A. asperula ssp. capricornu (Figs. 1, 6). Late season densities were several times lower and more broadly dispersed, with higher larval density areas in West Texas generally dominated by A. latifolia host records (Figs. 1, 6B). Early season monarch larval densities were lowest in parts of South and West Texas, while late season larval densities were lowest in parts of South Texas and Louisiana. The sample bias towards urban centers in the larval data supports our approach of focusing on the distribution of milkweed resources in the South Central US within the core habitats of the important focal milkweed species, for which occurrence data appeared to be less concentrated around urban centers (Fig. 7).

Milkweed niche modeling
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 confluence of the Pecos River and Rio Grande in Texas (Fig. 7B). Our data support 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).
Our models confirmed 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. 7A). 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 efforts to include more of a representative number of A. viridis occurrence data to the West and South (for discussion on niche model performance and feature selection, see Supplement, Discussion, Milkweed niche models).

Milkweed land cover selectivity
Negative land cover selectivity was found with Cultivated Crops for all four milkweed species. Waterbury et al. (2019) similarly found cultivated croplands had low utilization for northwestern milkweeds. We found a few A. oenotheroides in Cultivated Crops. Asclepias latifolia occurred in cotton fields of the Texas panhandle in the 1960's to 1970's before the increased use of glyphosate resistant cotton (S. Feagley, Soil and Crop Sciences Department, TAMU, retired, personal communication). The greater vegetative root budding of A. oenotheroides and the clonal nature of A. latifolia (Pellissier et al. 2016) probably make these species relatively more tolerant of cultivation. Widespread herbicide use could be lowering occurrence of A. oenotheroides and A. latifolia in Cultivated Crops in the South Central US, similar to the previously discussed reduction of rhizomatous A. syriaca in Midwest croplands (Zaya et al. 2017).
All four focal milkweeds had negative selectivity for Shrub/Scrub land cover, which could be related to lower soil moisture and anthropogenically lowered disturbance regimes fostering 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. Grassland Herbaceous land cover was of positive selectivity for A. viridis and A. oenotheroides, random selectivity for A. a. ssp. capricornu, and negative selectivity for A. latifolia (Fig. 8C). Waterbury et al. (2019) found highest utilization of northwestern milkweeds for grassland/herbaceous cover. Lower disturbance leading to less open areas in contemporary anthropogenically managed Grassland Herbaceous land covers may disfavor A. latifolia. Management practices with the potential to increase milkweed and nectar plant populations in grasslands include brush management, prescribed burning, prescribed grazing, and mowing (USDA NRCS 2015a, 2016a. The Developed-Low Intensity and Developed-Open Space land covers, which commonly included roadsides and urban areas, were of positive selectivity for two to three of the four milkweed species, respectively, except for A. latifolia. We found from 6 to 31% of our milkweed occurrences were associated with roadsides, depending on the species. 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). Management of roadside habitats is receiving increased focus for benefitting 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; Monarch CCAA/CCA Development Team 2020). Roadside habitats can be a threat to monarchs from roadkill (Kantola et al. 2019a;Mora Alvarez et al. 2019). Phillips et al. (2019) found roadside pollinator populations were higher a few meters away from the roadway, and Keilsohn et al. (2018) suggested that placing nectar plants a few meters from the roadway edge may reduce roadkill. Enhancement of milkweeds in urban gardens to benefit monarchs is also receiving attention (Fitzgerald 2015;Baker and Potter 2019). However, predation in 1 3 urban gardens can sometimes be higher than in rural habitats, potentially making them an ecological trap (Baker and Potter 2020; see also McGruddy et al. 2021).
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 inherent bias of iNaturalist observations along roadsides (Figs. S26-S29). However, there is still the potential for substantial error in the milkweed land cover selectivity analysis. For example, our sample bias reduction method could be overcompensating for bias of community science records near roads if a milkweed species had greater selectivity for roadside disturbed habitats than the larger target group of nectar plants (see also Ranc et al. 2017).

Seasonal milkweed resource areas
Previously discussed community science data biases probably reduce the accuracy of our milkweed resource area projections, and more studies with random sampling protocols are needed to refine predictions. However, our innovative combination of milkweed niche model core habitats with indices of seasonal landscape host selectivity and milkweed land cover selectivity provides valuable estimates of seasonal milkweed resource areas. The milkweed resource areas resolve major corridors providing monarch adult nectar and larval host plants for both the spring migration as well as the last summer generation and fall migration. The identified early and late season corridors of milkweed resources are similar and mostly encompassed within the previously identified Central and Coastal monarch fall migratory funnels of Tracy et al. (2019), except for Northeast Texas. The spring monarch migration is dispersed, and spring migratory routes are more difficult to identify than fall migration routes (Solensky 2004), which have been clearly defined (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 first 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, and most of the seasonal milkweed resource areas (Fig. 9). We identify a potential critical connectivity gap in milkweed land cover benefit resource areas in A. oenotheroides core habitat in South Texas (Fig. 9). Providing additional A. oenotheroides milkweed habitat in South Texas may improve spring and fall connectivity for monarch migration.
Disturbance in anthropogenic areas from mowing and burning favors growth of Midwest milkweeds, such as A. syriaca (Kaul and Wilsey 2019), as well as A. viridis (Baum and Sharber 2012;Baum and Mueller 2015) and probably other South Central US milkweeds. For example, regular late season mowing of grasslands in Oklahoma can increase the cover of A. viridis (Dee et al. 2016;Dee and Baum 2019). Leone et al. (2019) reported that monarchs were more abundant at sites managed by prescribed fire compared to grazing. Periodic mowing or burning is probably a critical disturbance factor needed to maintain open habitat for continued milkweed colonization and persistence in roadside and grassland habitats (Kaul and Wilsey 2019). Mowing of milkweed habitats at the wrong time during the spring and late summer/fall breeding seasons can reduce habitat availability for pollinators (Hopwood et al. 2015;Phillips et al. 2019;Waterbury et al. 2019), including immature and adult monarchs (Cariveau et al. 2020; Monarch Joint Venture 2020), and potentially reduce milkweed seeding and spread. Timing of mowing or burning in fields, roadsides, and urban open spaces can be especially important to maximize seasonal milkweed availability to 1 3 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.

Conclusion
Combined analyses of monarch community science data on landscape resource selectivity, distribution, and phenology can contribute to regional conservation assessments. Seasonal milkweed resource area estimates revealed variations in larval habitat quality within otherwise edapho-topo-climatically suitable milkweed habitats. Milkweed conservation in South Central US Grassland Herbaceous and Developed land covers, including roadsides and urban areas, has potential to benefit migrating monarchs by conserving or increasing critical milkweed nectar and host plant resources. Seasonal milkweed resource areas resolved important corridors in the Central and Coastal monarch migratory funnels where milkweeds provide adult nectar and larval food resources critical to functional connectivity for the spring and fall migrations through the South Central US. A potential critical gap in monarch migration functional connectivity was identified in South Texas. Seasonal milkweed resource area corridors and gaps can guide mitigation measures to conserve or restore monarch migratory habitats.

3
Ethical approval Not applicable.

Consent to participate Not applicable.
Consent for publication Not applicable.