Historic and Contemporary Land Use Shape Plant-Pollinator Networks and Community Composition

Globally, grasslands represent a critical but shrinking habitat for native plants and pollinators, with declines driven by alterations to landscape-scale habitat cover and local-scale disturbance regimes, among other factors. Specifically, as cities expand in size, an increasing proportion of regional pasture and grassland habitat is being replaced by urban development, and fewer periodic grazing and burning regimes are being supported locally, despite evidence that such regimes promote plant species richness and facilitate their interaction with native pollinators. The quantification of these plant-pollinator networks—through indices such as network connectance, specialization, nestedness, and robustness—can provide a unique opportunity to characterize key structural properties of species interactions and their response to human management and seasonal phenology. While urbanization and local disturbance regimes likely influence plant and pollinator communities and their interactions, past research in this area has primarily been conducted at limited spatial and temporal scales and has not typically quantified the impacts of both local and landscape forces on network properties. In this study, we investigate the effects of contemporary (past 10 years) and historic (prior 90 years) disturbance regimes on plant-pollinator community composition and network structure across more than 200 km of grassland in Central Texas. Our analyses indicate that for plant and pollinator communities, both contemporary and historic land management practices have led to significantly dissimilar community composition. Plant and pollinator richness and network nestedness are negatively correlated with phenological period, while pollinator richness is positively correlated with landscape-scale (2 km) urbanized land cover and is higher in historically grazed land, likely due to greater food and nesting resource availability. In contrast, we show that network connectance is positively correlated with phenological period and negatively correlated with landscape-scale urban cover. Finally, we show that pollinator robustness, a measure of resilience to plant species loss, is positively correlated with landscape-scale urbanization, likely due to greater redundancy provided by common weedy plant species. Overall, our results demonstrate that historic grazing regimes, current urbanization levels, and distinct phenological periods can simultaneously drive plant-pollinator community composition and network dynamics in shrinking but critical grassland ecosystems.


INTRODUCTION
In ecology, disturbances are often defined as stochastic processes which remove a substantial amount of biomass from an ecosystem at a given point in time (White and Pickett, 1985) and include processes such as wildfires, large mammal grazing, windstorms, and riverbank erosion . These disturbance regimes have dramatic impacts on vegetation density, ground cover, and soil substrate availability (reviewed in Mori, 2011), and can alter the diversity and flowering patterns of terrestrial plant communities (Sprugel, 1991;Moranz et al., 2012). Given that disturbance regimes can influence plant diversity and flowering patterns (e.g., Collins, 1987;Grundel et al., 2010;Moranz et al., 2012), they also have the potential to critically impact native insect pollinator communities [Peralta et al., 2017;reviewed in Koltz et al. (2018)], as many pollinators depend on plants for pollen and nectar to fuel their activity and provision thei1r brood (Michener, 2000). Specifically, contemporary land management practices such as prescribed burning, grazing, and mowing, which are employed globally (Bond and Keeley, 2005;Mapiye et al., 2008), have been shown to alter plant community composition (Collins, 1987;Howe, 1994) and ground cover (Gibson, 1988;Fidelis et al., 2012) across grassland systems. Indeed, there is growing evidence that current practices of low-level grazing (Vanbergen et al., 2014), occasional mowing (Weiner et al., 2011(Weiner et al., , 2014, and intermittent burning (Brown et al., 2017), may increase plant flowering, which could potentially alter pollinator abundance, diversity, and plant-pollinator interactions (e.g., Vanbergen et al., 2014). However, these studies are limited and have typically not considered the critical impact of historical land use when interpreting contemporary disturbance effects on plant and pollinator communities.
An increasingly large body of ecological research has demonstrated that historic land-use can have lasting impacts on contemporary vegetation assemblies (Greenlee and Langenheim, 1990;Floyd et al., 2003;Taverna et al., 2005;Ellis and Coppins, 2007;du Toit et al., 2016) and can even predict contemporary pollinator abundances (Cusser et al., 2015). For example, Taverna et al. (2005) found that current day vegetation patterns within hardwood tree stands were linked to past agricultural use, and Johnson et al. (2015) found that historic land-use, in the form of garden or building sites, differentially drove contemporary plant beta diversity. In a study of early successional forests, abandoned pastures supported very distinct contemporary plant species relative to abandoned crop fields (Benjamin et al., 2005). Likewise, in grassland-dominated bioregions, historic land use, including grazing or farming, has been documented as a significant predictor of pollinator abundances ∼50 years later (Foster et al., 2003), likely due to the fact that overwintering insects, like solitary bees, can have delayed population responses to pulsing floral resources that can be seen for many years after the disturbance event (Cusser et al., 2015). Within woody grasslands, past work has documented multi-decade lags in the impact of urbanization, where dissimilarities in contemporary vegetation composition were best predicted by landscape features ∼20-40 years before the survey period (e.g., road network density, percentage of natural area) (du Toit et al., 2016). In the same study, vegetation species richness of open grasslands was also predicted by recent (landscape features from 1 to 2 years before the survey) urbanization events. These studies have highlighted the importance of evaluating both contemporary and historic land management practices when quantifying drivers of plant and pollinator community composition and species interactions.
Species interaction networks, such as plant-pollinator networks, are excellent tools for quantifying the structure of mutualistic interactions [reviewed in Bascompte et al. (2003), Thébault and Fontaine (2010), and Dehling (2018)] and can also capture the impacts of local management and landscape composition on these critical interactions (e.g., Memmott et al., 2004). For example, networks can characterize the degree of connectedness (proportion of actual links to all possible links, sensu Dunne et al., 2002), specialization (the degree of niche partitioning across species, sensu Blüthgen et al., 2006), and nestedness (the degree to which specialized interactions are bound within more generalized interactions, sensu Bascompte et al., 2003) within a community. Interestingly, past studies have revealed that increases in floral species richness due to grazing can cause decreases in plant-pollinator network nestedness and increases in connectance (Vanbergen et al., 2014), indicating that communities with high biomass turnover rates may have a low buffer against specialized species loss. Indeed regional urbanization has been shown to lead to bird species loss, leaving only those species adapted to urban environments (Schneiberg et al., 2020); this may lead to increases in interaction evenness, where most animals are visiting the same plant species, typically because the remaining animal species tend to be generalist (Schneiberg et al., 2020). In addition, past work has indicated that decreases in local vegetation richness and structure can cause species loss and decreases in nestedness via a decline in floral resource availability (Moreira et al., 2015), though it is not known if this pattern persists across phenological periods within human-altered landscapes.
Indeed, because flowering duration for many plant species is short and pollinators are also often short-lived, the phenological stage of a plant and pollinator community is necessary to consider when quantifying interactions (Olesen et al., 2008), as network features such as network specialization can vary with the seasons and across different levels of floral resource availability (Harrison et al., 2018;Souza et al., 2018;Escobedo-Kenefic et al., 2020). Studies that have measured week-to-week variation in plantpollinator networks demonstrate that interactions can be highly variable in their connectance and robustness, leading to flexibility (changing values between timepoints) in the network structure relative to the cumulative network of all interactions in a season (CaraDonna and Waser, 2020). Connectance is another feature of networks that can vary based on the season, often due to fluctuations in the size of the network; networks tend to have lower connectance when there are more overall interactions occurring (Basilio et al., 2006). On the other hand, nestedness tends to increase with network size (Bascompte et al., 2003) and can be indicative of reduced interspecific competition and increased species coexistence (Bastolla et al., 2009). Despite the importance of phenology in network structure, past work conducted within human-altered landscapes has rarely explored the impacts of phenology alongside contemporary and historic land management practices.
In this study, we use a model grassland system in Central Texas to quantify the impacts of historic and contemporary land management on plant and pollinator community composition and interaction network structure. We hypothesize that similarity in contemporary disturbance regimes drives plant and pollinator community similarity more than historic land use. Specifically, we predict that sites that are currently actively managed (e.g., frequent fire or mowing) will exhibit greater plant and pollinator species richness than passively managed sites and will be more similar in composition than those with the same historic land use. Pollinators, and bees in particular, respond quickly to changes in vegetation and ground cover (Kimoto et al., 2012) thus, we anticipate that areas that are currently actively managed will have more similar plant and pollinator composition. We also hypothesize that landscapes with higher levels of surrounding urban cover will have lower network nestedness, robustness, specialization, and plant and pollinator species richness, due to a lack of colonizers that would otherwise contribute to community stability.

Study System
Research was conducted within 10 grassland study regions within rapidly urbanizing Central Texas, extending across more than 500 km (SW corner: 30.1015 N, 97.9608 W, NE corner: 33.6494, 95.6987 W) and comprising of state, federal, and nonprofit managed grasslands (described in Ritchie et al., 2016). To characterize each study region, five equidistant 50 m × 50 m plots were equally spaced along a 1.2 km linear transect within each region . Based on land manager interviews, we classified study regions based on both historic and current land management practices (similar to Gustavsson et al., 2007). Specifically, study regions were classified as either "historically farmed, " with hayfield or homestead farming between the late 1900's to the early 1990's, or "historically grazed, " with cattle and/or goats in the same time period. Contemporary management was categorized as either "active" or "passive, " with active being mowed or burned at least twice in the 10 years prior to sampling, and passive being unburned and unmowed in the 10 years prior to sampling (i.e., unmanaged). Management practices were orthogonal, as each current management category consisted of both historic land use types. Many past studies have grouped mowing and burning together given that they both remove substantial biomass and can have similar effects on plant production (e.g., MacDougall and Turkington, 2007;Dickson, 2019;Vermeire et al., 2020), though it is possible they distinctly impact individual pollinator species, a topic outside the scope of this study. We used the 2012 National Land Cover Database (30 m resolution, Homer et al., 2015) to characterize percent cover at a 2 km radius from the centroid of each study region, by first creating two broad categories, urban (made up of open, low, medium, and high development) (defined at mrlc.gov) and natural habitat (made up of grassland, forest, and shrub) (as per Plascencia and Philpott, 2017); because these categories were highly correlated (Pearson correlation = −0.878) we chose to use urban land cover in further analyses, as per many other studies (e.g., Matteson et al., 2013;Cusser et al., 2015;Plascencia and Philpott, 2017;Sexton and Emery, 2020). The 2 km radius was chosen to include typical pollinator foraging distances (Greenleaf et al., 2007), as in past studies conducted within the region (e.g., Ballare et al., 2019). In our study system, urban land cover was primarily comprised of low-level development in rural areas, including low-density housing in previously natural or agricultural areas (Hansen et al., 2005). This type of development is occasionally referred to as exurban and is one of the fastest types of land conversion in the United States, given substantial human relocation from cities to areas beyond the suburbs (Hansen et al., 2005).

Pollinator and Floral Resource Sampling
In the summers of 2012 and 2013, we surveyed floral communities via quadrat surveys, and we sampled native bees and butterflies via netting and trapping at each plot during three distinct phenological periods: early bloom (April 18 to May 15), mid-bloom (May 20 to June 16), and late bloom (June 20 to July 14) (Ritchie et al., 2016). Plots were sampled once during each period. Specifically, during each of the three phenological periods, we measured floral species richness and floral density in 30 1 m × 1 m quadrats per site. The quadrats were evenly positioned 4 m apart along three 50 m transects running from North to South that were located at 10 m, 25 m, and 40 m from the NW corner of the plot (Ritchie et al., 2016). Specifically, within each quadrat, the number of forb inflorescences per species were counted. We also quantified ground cover (bare ground, vegetation, rocky, and impervious cover) within each quadrat, and across the study system, and we measured the size of five flowering heads per species to calculate total floral cover. Ground cover metrics were highly correlated (e.g., bare ground and vegetation cor = −0.434, p < 0.0001), therefore we focused on bare ground for further analyses (Ballare et al., 2019).
For the pollinator surveys, 2 researchers netted for 30 min by walking slowly back and forth on the east or west side of the plot for 15 min and then switching sides with their partner for the other 15 min, between the hours of 7 am and 12 pm during each phenological period, and only on sunny days. Timers were not stopped while insects were put into kill vials (Ballare et al., 2019). During this 30-min time period, all native bees and butterflies observed foraging on flowers were caught and the flower species was recorded, as in past plant-pollinator network studies (Winfree et al., 2014). We focused on bees and butterflies as in previous studies (Buhk et al., 2018;Librán-Embid et al., 2021) and because these groups are among the most common and effective pollinator taxa in the study region (Sexton and Emery, 2020). Individuals were placed in separate kill jars and all individuals were pinned, labeled, and identified to species. After completing visitation surveys, pan trapping was conducted by placing 30 pan traps (6-oz plastic bowls, SOLO model number: PB9-0099) 1 meter apart and alternating by color between white, blue, and yellow, in an X-formation in the middle of the study plot (LeBuhn et al., 2003) and blue vane trapping was conducted by placing 4 blue vane traps in the center of the plot, hanging one meter off the ground on a wooden stand (Ballare et al., 2019). The pan-traps were filled with 4 oz of a diluted soap water (1 gallon water: 1 tbsp Dawn dishwashing soap) and left in the field for 24 h, after which insects were stored in 90% ethanol, and blue vane traps were left in the field for 5 days before specimens were collected and stored in 90% ethanol (as per Ballare et al., 2019). To comprehensively characterize each study region, the ground cover data, netted pollinator and plant data, and trapped pollinator data from the five plots within each region were combined (Baldock et al., 2015) for each of the phenological periods (similar to Prendergast and Ollerton, 2021), for a total of 3 ground cover, netted pollinator and plant, and trapped pollinator datasets per study region per year (n = 60 per dataset).

Network Analysis
Plant-pollinator networks were created from the netted data using the R package bipartite (Dormann et al., 2009). We focused on network-level nestedness (NODF), specialization (H2), connectance, higher (pollinator) and lower (plant) level species richness, and higher-and lower-level robustness. Nestedness is a term used to describe the structure or organization of network interactions, where more generalist species from both the higher and lower orders (animals and plants) interact with more specialist species (Bascompte et al., 2003) and is measured by the overlap and decreasing fill of the plant-pollinator matrix, with values between 0 and 100 where higher values indicate higher nestedness (Almeida-Neto et al., 2008;Ulrich et al., 2009). Past work suggests that greater nestedness is indicative of lower interspecific competition and greater coexistence between species (Bastolla et al., 2009) and thus is a signature of more stable communities (sensu May, 1972) that rebound more quickly to equilibrium following perturbations (Thébault and Fontaine, 2010). Specialization (H2) describes the entire network's level of specialization, or the degree of niche partitioning across species (Blüthgen et al., 2006) compared to the expected interactions given the number of interacting species (Dormann et al., 2009), where a value of 0 represents no specialization and a value of 1 represents a completely specialized network. Connectance is the proportion of all possible links in the network that are connected, where a value of 0 indicates no interactions and a value of 1 indicates that all plant species are interacting with all pollinator species. Higher levels of connectance lead to higher levels of robustness to extinction (Dunne et al., 2002) as well as greater stability (sensu May, 1972). Higher level robustness characterizes the pollinator guild and lower level robustness characterizes the plant guild, where the area below the extinction curve quantifies the robustness of the system to species loss; this is based on the assumption that if a fraction of species from one guild are eliminated, then many species of the other guild will go extinct (Dormann et al., 2008). For nestedness, robustness, and specialization indices, we converted these values into Z-scores by calculating the mean of the network index divided by the standard deviation of 1,000 null models created with the nullmodel function in the package bipartite and subtracting this from the observed network level value (Vázquez and Aizen, 2003). We decided to use the z-scores for this because it allows comparison of the focal network to what is expected if all interactions were random (Gotelli, 2001).

Composition Analyses
We tested for differences in plant and pollinator community composition across historic land-use, current management types, and phenological period using permutational MANOVA using the adonis function in the R package vegan, with year as a controlled stratification factor given that the same regions were sampled 2 years in a row. Specifically, PERMANOVAs were conducted on the raw, log(X + 1) transformed, and presence absence data to control for compositional differences driven primarily by abundance (Ballare et al., 2019). We used bray-curtis dissimilarity as this metric is commonly used for community studies (Burkle and Alarcón, 2011). We also ran the PERMANOVAs using the morisita-horn metric and found similar results (Supplementary Table 5). We used non-metric multidimensional scaling with the metaMDS function in the R package vegan to visualize differences in the communities.

Habitat Indicator Species Analysis
We used the multipatt function in R to perform multi-level pattern analysis in the indicspecies package to quantify indicator species for both contemporary management and historic land use types (package "indicspecies"). Indicator species capture the strength of the relationship between species and the groups of regions where they occur, and indicate which species are the predominant species in that habitat type and not in others (Cáceres and Legendre, 2009).

Regression Analyses
We used regression models to investigate the impact of five predictor variables: historical land use (farmed or grazed), current land management type (active or passive), landscapelevel urban habitat cover (2 km radius), phenological period (1, 2, or 3), and local bare ground cover on two response variables from the trapped datasets, pollinator abundance and richness, and seven network response variables from the netted data: higherlevel species richness, lower-level species richness, specialization, connectance, nestedness, higher-level robustness, and lower-level robustness. All continuous predictor variables (bare ground and percent grassland habitat cover) were scaled, and the year and study region were used as random effects in all models. We tested for an interaction effect between urban cover and phenological period and found no significant effect in any models, and therefore decided not to include the interaction in the final models. We created generalized linear mixed effect models using the glmer function in the R package lme4 with a Poisson distribution for higher-and lower-level species richness, given these are count data. We created generalized linear mixed models using the glmmTMB function in the R package glmmTMB for specialization and connectance with a beta distribution, because they range from 0 to 1. Nestedness and higher and lower-level robustness were normally distributed so we used linear mixed effects models using the lmer function in the R package lme4.
All models were checked for collinearity by calculating a variance inflation factor (VIF) using the car package in R (Fox and Weisberg, 2019) and all models were below our conservative cut-off of 3. Finally, AICc-based model selection was run for all models using the dredge function in the R package MuMIn, given that AICc is particularly suitable for smaller datasets (Bedrick and Tsai, 1994), and we used a delta AICc of 2 for averaging top models within this bound.

RESULTS
We recorded and identified a total of 223,632 inflorescences in the vegetation surveys, with between 1 and 16 plant species represented in each study region per phenological period (mean 7.82 SE 0.45). Bare ground covered 7.38% of the average surface in each study region. We collected 16,950 insects in the pan and blue vane traps, consisting of 240 different pollinator species, ranging from 10 to 67 species in each study region per phenological period. The three most abundant bee species found in the pan-and blue vane-traps were Lasioglossum TX. sp.3 (2,839 individuals), Lasioglossum coactum (1,674 individuals), and Lasioglossum bardum (1,587 individuals); the three most abundant butterfly species were Lerodea eufala (284 individuals), Pyrisitia lisa (185 individuals), and Pyrgus communis albescens (101 individuals).
We observed a total of 2,655 total interactions in the netted surveys, ranging from 1 to 16 plant species and 2-32 pollinator species in each study region per phenological period (mean 11.58, SE 0.43), and a total of 177 pollinator species (bees and butterflies) and 112 plant species overall. The most abundant bee species found in the netted surveys was Bombus pensylvanicus (244 interactions), followed by Xylocopa virginica (223 interactions) and Melissodes coreopsis (141 interactions). The most abundant butterflies caught were Euristrymon Ontario (97 interactions), Euptoieta Claudia (57 interactions), and Nathalis iole (54 interactions).

Composition Analyses
Our PERMANOVAs showed that the floral communities were significantly different between historic land-use types for the raw (P = 0.001), log-transformed (P = 0.001), and presence absence data (P = 0.001) (raw data visualized, Figure 4A) and significantly different between current management groups for the raw (P = 0.015), log-transformed (P = 0.023), but not the presence absence data (P = 0.089) (raw data visualized, Figure 4B). Results were nearly identical when using the Morisita-horn method in place of the Bray-Curtis method (Supplementary Table 4). Trapped pollinators were also significantly different between historic land-use types for the raw (P = 0.001), log-transformed (P = 0.001), and presenceabsence data (P = 0.001) (raw data visualized, Figure 4C) and significantly different between current management groups for the raw (P = 0.002), log-transformed (P = 0.001), and presenceabsence data (P = 0.034) (raw data visualized, Figure 4D). Netted pollinator communities were significantly different between historic land-use types for raw (P = 0.001), logtransformed (P = 0.001), and presence-absence data (P = 0.001) (raw data visualized, Figure 4E), while current management was not a significant driver of different netted pollinator communities with the raw (P = 0.078), log-transformed (P = 0.164), or presence absence data (P = 0.251) (raw data visualized, Figure 4F). For both communities, netted and trapped pollinators, results were identical when using the Morisita-horn method in place of the Bray-Curtis method (Supplementary Table 4).

Habitat Indicator Species
Within the trapped pollinators, out of a total of 235 species, 13 were found significantly more in the historically farmed sites, 22 species were found significantly more often in the historically grazed sites, 4 species were found significantly more often in contemporary actively managed sites, and only 1 species was found more often in contemporary unmanaged sites. Of these, the top five indicator species found significantly more frequently in historically farmed sites were Lasioglossum disparile Within the netted pollinators, out of a total of 177 species collected, 5 species were found significantly more in historically farmed sites, 8 were found significantly more in historically grazed sites, and 3 were found more in actively managed sites than passively managed sites, which had no prevalent indicator species. The indicator species found significantly more Hylephila phyleus (stat = 0.408, p = 0.010), and Bombus fraternus (stat = 0.408, p = 0.038), while the top five species found significantly more frequently in historically grazed sites were Megachile policaris (stat = 0.735, p = 0.001), Pyrgus communis albenscens (stat = 0.693, p = 0.002), Diadasia rinconis (stat = 0.577, p = 0.004), Nathalis iole (stat = 0.576, p = 0.016), and Lasioglossum coactum. The three species found significantly more in contemporary actively managed sites were Lasioglossum disparile (stat = 0.613, p = 0.034), Junonia coenia (stat = 0.548, p = 0.006), and Bobmus griseocollis (stat = 0.429, p = 0.027) (Supplementary Table 5).

DISCUSSION
In this study, we found that pollinator richness and abundance were significantly higher in historically grazed vs farmed land and that both floral and pollinator richness decreased with phenological period. In addition, we found that while network connectance was higher in habitats with contemporary active management, robustness was lower, indicating that these mowing and burning management practices may increase species interactions but may not necessarily strengthen plantpollinator network stability. Pollinator richness and robustness were also higher in landscapes with higher surrounding urban cover, indicating that moderately developed spaces may provide novel resources for pollinators. Finally, we found striking compositional differences in floral and pollinator communities based on both contemporary and historic land use practices.

Plant and Pollinator Richness
We found that both floral species richness and pollinator species richness was significantly greater in grassland areas that were historically grazed. Field experiments conducted within grasslands have shown that grazing, and the addition of burning to grazed plots, can increase plant species richness [Collins, 1987;Gibson, 1988;reviewed in Valkó et al. (2014)]. This is because grazing is a gradual disturbance process that removes biomass from grassland systems, allowing for colonists and seeds within the seed bank to establish. Disturbances increase environmental heterogeneity by changing soil characteristics, like nitrogen availability (Baer et al., 2016), allowing species coexistence within different patches , resulting in patches of grassland that contain unique species composition (e.g., Collins, 1987). This benefit to the flowering plant community can have cascading effects on arthropod diversity (van Klink et al., 2015). A recent meta-analysis summarizes substantial past work that resonates with our findings, where Hymenoptera respond positively to wildfire disturbance (Carbone et al., 2019), indicating an important role of biomass removal for bee and wasp pollinators.
In our study, we also found greater pollinator richness in grasslands with greater urban land cover, perhaps driven by greater floral and nesting resource availability in these landscapes (reviewed in Sexton and Emery, 2020). Indeed, a variety of nesting and food resources are necessary for many pollinators who forage in open grasslands but nest FIGURE 4 | NMDS plots showing the differences between communities for plants between the historic land-use types (A) and current management types (B), for trapped pollinator communities between historic (C) and current (D), and for netted pollinators for historic (E) and current practices (F).
in wood or leaf litter (Cane and Tepedino, 2001), and past studies have documented cases where urban landscape can provide some of these critical resources (Cane et al., 2007). Even within highly urbanized areas, studies have shown that low levels of developed land can increase floral diversity and nesting habitat at a landscape-scale (Matteson et al., 2013), potentially allowing for a greater diversity of pollinators to persist in neighboring natural habitat. Lowenstein et al. (2014) found that more dense neighborhoods supported a greater diversity of flowering plants, leading to greater bee diversity in sites with greater human population density. Residential and community gardens have also been shown to be pollinator abundance hotspots due to high floral resource availability (Baldock et al., 2019). Finally, more urbanized landscapes often exhibit higher levels of habitat heterogeneity (McDonnell and Pickett, 1990) and this heterogeneity may lead to increases in diversity in many arthropods (Báldi, 2008). In our study system, urban cover was comprised primarily of low and open development (i.e., large-lot single-family homes, parks, golf courses) and this type of land use is increasing across the southern United States at rapid rate (Johnson and Beale, 1998). Outside of exurban systems, within agroecosystems, a recent meta-analysis that modeled the effects of land-use intensity on pollinator biodiversity also found that low levels of land-use intensity (e.g., minimal-use urban) are indeed beneficial for pollinator biodiversity (Millard et al., 2021), suggesting an important role of promoting habitat heterogeneity in ever-expanding lowdevelopment exurban landscapes (Cusser et al., 2018). We also found that pollinator and floral species richness decreased significantly throughout the flowering season, from May to July, as the temperature in the Southern Plains tend to peak (Griffiths and Strauss, 1985) and the abundance of wildflowers tend to decrease. This pattern has been documented in many temperate grassland systems (Wilsey et al., 2011) as well as in our central Texas study system (Ritchie et al., 2016;Ballare et al., 2019). Such pronounced shifts in floral and pollinator richness likely influence pollinator foraging preference, as previously documented in in our study system (Ritchie et al., 2016), but also impact plant-pollinator network structure.

Network Characteristics
Indeed, we found that several network indices, including connectance, nestedness, and robustness were significantly negatively correlated with phenological period. While some studies have found similar patterns with connectance (Vázquez et al., 2009;CaraDonna and Waser, 2020) others have found the opposite relationship (Basilio et al., 2006). The studies that show similar patterns, where connectances decreases as the number of plant and pollinator species decreases, all sampled throughout the flowering season to account for variability in network structure. Specifically, CaraDonna and Waser (2020) describe interaction flexibility between species, where network structure changes from week to week, resulting in differences in connectance, nestedness, and specialization, over the summer growing season. CaraDonna and Waser (2020) also found that throughout the phenological season (where floral species richness declined), both connectance and nestedness were dependent on species richness. Interestingly, we found that nestedness (NODF) was negatively correlated with phenological period, meaning that the networks became less nested over the flowering season. This same pattern has been documented in other grassland system (Dunne et al., 2002), suggesting that networks are more robust when there are more species present. In past simulation studies, robustness has indeed been positively correlated with topological plasticity, indicating that when species are removed from the system, the network is more robust if the species are able to fill in the roles of the newly extinct species (Somaye et al., 2020).
We also found that increases in urban cover at a landscapelevel led to increases in higher level robustness; meanwhile, connectance and higher-level robustness were lower and higher in passively managed regions, respectively. Our finding of positive relationships between robustness and urban cover could indicate that developed spaces can increase beneficial habitat and support disparate species, as seen in pollinator studies conducted in more urbanized landscapes (Matteson et al., 2013). While some network studies conducted in wetter regions have not documented this pattern (Moreira et al., 2015), our focus on more water-limited ecoregions, highlights the potential for urban habitat to provide floral resources not found in natural habitat within arid landscapes (Cane et al., 2006;Baldock et al., 2015;Wenzel et al., 2020). We also found that connectance was lower with contemporary passive land management and decreased with increasing regional exurban cover. Because connectance captures the proportion of realized links and can predict network stability (Poisot and Gravel, 2014), our finding indicates land management techniques such as grazing can impact stability, as seen in past studies (Vanbergen et al., 2014). Overall, our findings demonstrate that network structure between plant-pollinator interactions is primarily a function of contemporary disturbances.

Plant and Pollinator Community Composition
We found that plant and pollinator community composition were driven by both historic and contemporary land use practices, as seen in a number of past studies on plants and pollinators (Benjamin et al., 2005;Cusser et al., 2018). Indeed, focusing just on grassland systems, both contemporary (Carvell, 2002) and historic land-use (Cusser et al., 2018) also played a major role in determining species richness and community composition (Howe, 1994;Foster et al., 2003;Ellis and Coppins, 2007;du Toit et al., 2016). In other words, longterm disturbance regimes, such as wildfire burns, can invoke differences in plant species assemblages, even in the presence of contemporary burn or mowing practices (Fidelis et al., 2012), and these "land-use legacies" may be quite strong, though they are often overlooked when focusing on contemporary species patterns (reviewed in Perring et al., 2016). For example, plant communities (Mattingly and Orrock, 2013), lichen (Berglund and Jonsson, 2005;Ellis and Coppins, 2007), and even other flying insect classes, such as hoverflies and butterflies (Sang et al., 2010;Bommarco et al., 2014) have exhibited significant compositional responses to historic land-use regimes. In another study, Sang et al. (2010) found that species richness of habitat specialist butterfly and moth species was positively related to both current and historic surrounding natural habitat area within a 2 km radius (Sang et al., 2010). Our study similarly indicates that both contemporary and historic management have strong impacts on plant and pollinator community composition and species interactions.

Habitat Indicator Species
We found that more than 40 species of pollinators were significantly driving differences in community composition between historically farmed and grazed sites and 8 species significantly driving differences between contemporary actively managed and passively managed sites, likely due to differences in the landscape's ability to provide resources for bees differing in body size and nesting preferences (Ballare et al., 2019). There were more pollinator species differentially representing the historically grazed sites, 26 unique species between netted and trapped methods combined, while the historically farmed sites had 16 unique species driving differences in community composition. The indicator species for the historically grazed sites tended to be more heterogenous in nesting type and foraging specialization (combination of oligo-and polylectic species) while indicator species in the historically farmed sites tended to be more polylectic. These differences in communities could be due to lasting effects that farming and large-mammal grazing have on the respective soil and vegetation types (Foster et al., 2003;Tappeiner et al., 2020). Specifically, grazing may be more similar to the natural history of disturbances in the region (Frank et al., 1998), where specialized oligolectic pollinator species may have adapted over time. Out of the five netted indicator species that differentiated farmed from grazed sites, three of them were from the genus Bombus. Widely known for their large body size, sociality, and surface or underground colony nests, these bumblebee species may be doing particularly well in the historically farmed sites because of the potential to use abandoned small mammal burrows as nesting resources (Kells and Goulson, 2003). Bumblebees have been recorded to use abandoned rodent holes as the base for their nest (Harder and Real, 1987), and rodents often live near human-dominated landscapes (Purvis et al., 2020). The greater representation of this group suggests that conversion of historically farmed land back to grassland habitat may be particularly beneficial conservation strategy for this group, which is believe to be declining more than other species.
The indicator species that distinguished the actively managed contemporary habitats from the passively managed ones were diverse in body size, with 4 out of the 7 species being Lasioglossum and one from the genus Bombus. The only significant indicator species found in the passively managed habitats was the butterfly Pyrgus communis albenescens, which was also significantly found in the historically grazed sites.

Conclusion
Overall, we found that land use history and contemporary land management can differentially impact the community composition and species interactions of plants and pollinators in grassland ecosystems. Specifically, we show that both historic (∼80 years prior) grazing and farming practices as well as contemporary (∼10 years prior) burning and mowing have lasting impacts on the composition of plant and pollinator communities. Our results also demonstrate that landscape-level urban cover is a driver of network structure and may lead to higher levels of plant-pollinator network robustness, especially in arid grassland systems. Additionally, our results show that plant and pollinator richness decrease across phenological periods, resulting in altered network structure. Finally, we show that indicator species which characterize historically grazed sites exhibit substantial diet and nesting heterogeneity, likely driven by similarity between pastoral systems and natural disturbance regimes in the bioregion.
Our study shows that to properly evaluate an area for conservation efforts, historic land use should be considered as this can have a lasting impact on current communities. Furthermore, future research should also not neglect this factor in community ecology work. Future research directions include investigating the impacts of other disturbance regimes, such as wildfires, and the impact that a combination of disturbances in one area has on the surrounding communities.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
CC conducted analysis and interpreted the data, made figures and tables, and wrote the manuscript. SJ designed the study and reviewed and revised the manuscript at all steps. JN was responsible for all taxonomic identification and read and suggested feedback on the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
CC was funded as a graduate student through the UT Austin Integrative Biology department, Mentoring Dean's Fellowship and the IB research grants.