Invertebrate and Plant Community Diversity of an Illinois Corn–Soybean Field with Integrated Shrub Willow Bioenergy Buffers

: Perennial bioenergy crop production within intensively managed agricultural landscapes has the potential to improve the sustainability, resiliency, and diversity of these landscapes. Perennial crops are ideal because of their high production potential on marginal lands relative to grain crops (e.g., corn and soybean) and their ability to provide additional ecosystem service beneﬁts. When agricultural landscapes are designed to target speciﬁc services, determining the non-targeted services of perennial bioenergy crops can further promote their adoption. This 3-year study addresses this proposition by evaluating the canopy invertebrates and understory plant (non-target crop) communities using bee bowls and point measurement of ground coverage, respectively, within a grain ﬁeld integrated with shrub willow buffer systems designed for nutrient loss reduction. Greater plant diversity and richness were observed under willow than under grain, resembling that of the surrounding riparian community with more perennial, native species. However, the same relationship did not hold true for invertebrates, with seasonality having a signiﬁcant inﬂuence resulting in similar communities observed in willow and grain plots. The presence of unique plant and invertebrate species in both willow and grain crops as well as foraging pollinators on both crop and non-target crop species highlights the importance of habitat heterogeneity for supporting biodiversity and the potential beneﬁts of buffer bioenergy landscape designs.


Introduction
Agricultural landscapes provide many human services, including the production of food, feed, and fiber, and, more recently, feedstocks for bioenergy and bioproducts. However, intensively managed modern agricultural systems negatively impact the natural ecosystems' support services, such as hydrological services, soil quality, and pollination. These production systems rely on support services for long-term productivity and resiliency [1][2][3]. Biodiversity, in particular, has been highly impacted by agricultural intensification, which has resulted in increased use of chemicals and habitat loss [2,4,5]. The loss of biodiversity can have significant economic impacts, with pollination services estimated to contribute USD 217 billion to the global economy [4,6] and natural pest control estimated to save U.S. agriculture USD 13.6 billion yr −1 [1,5].
The promotion of perennial bioenergy crop production (such as perennial grasses and woody crops) for biodiversity support is attractive, as perennial bioenergy crops generally require less pesticide and fertilizer usage than traditional crops, resulting in a lower environmental impact compared to grain crops such as corn (Zea mays L.) and soybean (Glycine Max [L.] Merr) [3,[7][8][9]. The presence of bioenergy crops within the agricultural

Experimental Design
The biodiversity study was conducted on a 6.5-ha corn-soybean field in Fairbury, Illinois, U.S.A. (40.74 • N, 88.50 • W) in 2017 to 2019 ( Figure 1). This study was a component of a larger landscape design study that initially focused on developing a nutrient loss reduction and recovery system. The site was established in 2013 and was described in Ssegane et al. [23] and Zumpf et al. [24]. The average annual precipitation for the site from 2010 to 2018 was 1454 mm (at the Fairbury wastewater treatment plant; [25]), and the long-term average for the area was around 1535 mm (from nearby Dwight, 1981-2010; [26]). The field site is neighbored by Indian Creek to the west, residential and commercial areas to the north and south, and corn-soybean production to the east ( Figure 1).
Three habitats, or landcovers, were monitored for this study, including grain, willow, and riparian edge habitats. For simplicity of terminology, they will be denoted as landcovers. The riparian edge habitat along Indian Creek was included as an unmanaged check, as it was expected to be impacted by neighboring production. However, there was no direct management of those areas. Riparian plots were classified as edge habitat, as traps and vegetation monitoring were only a few meters away from the production edge.
The field experiment was arranged as a randomized block design with two replications of grain (corn and soybean) and short-rotation shrub willow (Salix miyabeana Seemen 'SX61') in 0.2-ha plots on two different soil types (marginal and non-marginal land). The northern plots (N-plots) were placed on a lowland floodplain with Comfrey loam soil (fine-loamy, mixed, superactive, mesic Cumulic Endoaquolls; 0-2% slope), which is more fertile with a higher crop productivity for corn. The southern plots (S-plots) were placed Sustainability 2021, 13, 12280 3 of 19 on a marginal, uphill plain with Symerton silt loam (fine, loamy, mixed, superactive, mesic Oxyaquic Argiudolls; 2-10% slope), which has greater soil erosion and susceptibility to nutrient loss [24]. The field experiment was arranged as a randomized block design with two replications of grain (corn and soybean) and short-rotation shrub willow (Salix miyabeana Seemen 'SX61') in 0.2-ha plots on two different soil types (marginal and non-marginal land). The northern plots (N-plots) were placed on a lowland floodplain with Comfrey loam soil (fine-loamy, mixed, superactive, mesic Cumulic Endoaquolls; 0-2% slope), which is more fertile with a higher crop productivity for corn. The southern plots (S-plots) were placed on a marginal, uphill plain with Symerton silt loam (fine, loamy, mixed, superactive, mesic Oxyaquic Argiudolls; 2-10% slope), which has greater soil erosion and susceptibility to nutrient loss [24].

Planting and Management
Willow plots were planted at a density of 15,300 cuttings ha -1 in the spring of 2013 using a double-row planting design [27]. N-plots had a lower plant density due to equipment calibration issues, and hand plantings to fill gaps were unsuccessful. Willow plots were not directly fertilized, but instead, the leached nutrients from upslope grain served as their nutrient source. Weeds were managed by mowing after the establishment period due to foliar damage from the application of pendimethalin (5.7 L ha -1 ) and Poast (sethoxydim: 2.4 L ha -1 ) herbicides in the first year. Willow plots were then left unmanaged after 2017. Willows were coppiced after the first season of growth, and harvest occurred on a three-year cycle. Sampling in this study was conducted across the three years of the second harvest cycle.
Grain plots had a soybean (no-till)-corn (till)-soybean (till) rotation during the study period (2017-2019). Historically, the field was under continuous corn (till) from 2009 to 2016. Weeds during soybean years were managed with a pre-emergent herbicide (Boundary: S-metolachlor and metribuzin) and post-emergent herbicides (Flexstar: omesafen at 1.5 L ha -1 and Fusilade: fluazifop-p-butyl at 0.6 L ha -1 ). No nitrogen (N) fertilizer was applied during the soybean years, but corn was managed with 248 kg N ha -1 yr -1 . Weeds

Planting and Management
Willow plots were planted at a density of 15,300 cuttings ha −1 in the spring of 2013 using a double-row planting design [27]. N-plots had a lower plant density due to equipment calibration issues, and hand plantings to fill gaps were unsuccessful. Willow plots were not directly fertilized, but instead, the leached nutrients from upslope grain served as their nutrient source. Weeds were managed by mowing after the establishment period due to foliar damage from the application of pendimethalin (5.7 L ha −1 ) and Poast (sethoxydim: 2.4 L ha −1 ) herbicides in the first year. Willow plots were then left unmanaged after 2017. Willows were coppiced after the first season of growth, and harvest occurred on a three-year cycle. Sampling in this study was conducted across the three years of the second harvest cycle.
Grain plots had a soybean (no-till)-corn (till)-soybean (till) rotation during the study period (2017-2019). Historically, the field was under continuous corn (till) from 2009 to 2016. Weeds during soybean years were managed with a pre-emergent herbicide (Boundary: S-metolachlor and metribuzin) and post-emergent herbicides (Flexstar: omesafen at 1.5 L ha −1 and Fusilade: fluazifop-p-butyl at 0.6 L ha −1 ). No nitrogen (N) fertilizer was applied during the soybean years, but corn was managed with 248 kg N ha −1 yr −1 . Weeds during corn production were managed with Resicore (acetochlor, mesotrione, and clopyralid at a rate between 2.9 and 3.5 L ha −1 ) and Atrazine (applied between 1.2 and 3.5 L ha −1 ). In the case of any carryover effects from prior years, it should be noted that continuous corn were managed with pre-emergent (3.5 L ha −1 of Harness Xtra: acetochlor) and post-emergent (0.04 L ha −1 of Impact: topramezone methanone) herbicides.

Invertebrate Sampling
Canopy invertebrate samples were collected using colored bee bowls (19 cm diameter and 3.8 cm deep). Three colors (blue, yellow, and white) were initially used [28]; however, since similar invertebrate communities were found between white and yellow bowls in 2017, only white and blue bowls were used in 2018-2019 to increase sampling efficiency. The bowls were placed on a 1-m-high platform and were filled with a solution of deionized water and Dawn dish soap (which acted as a surfactant: Procter & Gamble Corp., Cincinnati, OH, USA). Each platform was placed at a minimum of 5 m away from each other in the plots to avoid interference [28], with three platforms placed in each plot. Invertebrate sampling targeted flowering periods of plant and crop species across the three landcovers (April, June, July, and September). The bowls were left in the field for 24 h on sunny, calm days. Collected invertebrates were preserved in a water and ethanol solution and were identified to taxa and functional group based on previous literature/knowledge of the species.

Understory Vegetation Sampling
Understory plant species' percent ground cover was measured in 2018 and 2019. A 1-m 2 PVC pipe quadrat was used to visually estimate the absolute cover of vegetation species present, excluding the main crop (in the case of willow and grain plots). In 2018, only one sampling area immediately below each set of bee bowl platforms was monitored (platform served as the middle of the 1 by 1-m plot). The sparse plant coverage during different times of the years, however, resulted in a greater number of quadrants with zero plant coverage, which underrepresent the existing plant community within the plots. As a result, in 2019, the number of understory plant sampling quadrats was increased to three around each bee bowl platform, each approximately 1.5 m out from the base of the platform in a triangular pattern. The resulting plant coverage was then averaged across the three quadrats per platform.

Diversity and Species Richness
Relative invertebrate abundance was calculated for each sampling location (colored bowls were combined) on each sampling date. The Shannon-Wiener Diversity Index was calculated for each sampling location and was averaged across landcover types. The influence of landcover type on invertebrate diversity and richness was analyzed using linear (LM) and generalized linear models (GLM), respectively, in R statistical software using the "vegan" package [29,30], with landcover type and sampling date serving as fixed factors. Poisson regression with log transformation was used in the GLMs for analyzing invertebrate richness. The "poly" function was used to determine if quadratic trends for sampling date were to be included in the models (p < 0.05). Analysis of variance (ANOVA; alpha = 0.05) was then used to identify significant model factors, with nonsignificant factors being dropped to simplify the models. Pairwise comparisons with a Tukey adjustment were used to test within-factor differences. Least-square (LS) means were back-transformed from log scale if GLMs were used.
For the understory plant community, the absolute cover of each plant species was converted to relative cover by dividing the percent cover of each species by the total plant coverage in the 1-m 2 quadrat [31]. Shannon-Wiener's Diversity Index was calculated using relative cover as a measure of abundance [32]. LMs and GLMs were also used to assess the effect of landcover and sampling date on plant community diversity and species richness. Following Fletcher et al. [33], two models were used to analyze plant diversity in 2018 because of the large number of 0s (plant absence) in the dataset. First, a GLM with the binomial family and logit link were used with a presence/absence model (1 s or 0 s) [33]. The second model was an LM, which used a dataset with the 0s removed. In 2019, only the LM model was needed due to the increased sampling. Significant LM and GLM fixed factors were determined using an ANOVA with within-factor differences evaluated with LSmeans pairwise comparisons. Plant species richness was analyzed the same way as invertebrate richness in both years.

Community Composition
A principal component analysis (PCA) was used to assess invertebrate and plant community composition. Any species observed only once was removed from the dataset, and each sampling year was run separately. Significant PCA axes were determined using the Kaiser's rule of 1. A k-means cluster analysis was used to identify groupings among sample communities using the R package "MASS" [34]. Scree plots were used to determine the number of clusters. The community composition associated with landcover and sampling date was further assessed with a random forest (RF) model ("randomForest" R package) using landcover and sampling date as the response variables and invertebrate species abundance as the predictor [35]. The "caret" and "e1071" packages were used to determine the optimal number of variables to try at each split [36,37], and the number of trees was determined using error plots [38].

Plant and Invertebrate Community Relationship
A symmetric and predictive co-correspondence analysis (Co-CA) was used to assess the relationship between the invertebrate and vegetation communities with the "cocorresp" R package [39]. Analysis followed Simpson [40] using raw insect species counts and vegetation percent coverage converted into 1-9 values on the van der Maarel scale according to van der Maarel [41]. Insect data were log-transformed using the "log1p" function. Matrix symmetry for the symmetric Co-CA was determined by how much of the total matrix variation was explained by the Co-CA axes and how the proportion of variation explained compared between the two matrices. For the predictive Co-CA (which used the vegetation matrix as the predictor matrix), a cross-validation fit percent method (using the "crossval" function; axes values > 0) and a permutations test (99 permutations; p-value > 0.5) were used to check if the model fit was significant. If the methods disagreed, then the lower number of significant axes of the two methods was taken [42].

Invertebrate Community
The sampling date was found to have a larger influence on the invertebrate community diversity and richness than the landcover type. During the soybean years (2017 and 2019), invertebrate diversity was not found to differ between landcover types (p = 0.09 and 0.34, respectively). However, the species diversity significantly increased over time across landcover types (Figure 2a,c) with a 10% increase from July to August in 2017 (p = 0.03) and a 20% increase in invertebrate diversity from April to July in 2019 (p < 0.03 for all comparisons). A general decline in species diversity was observed between July and September of 2019; however, that change was not significant (Figure 2c; p > 0.74). Species richness followed similar patterns, with a 25% increase across landcover types in 2017 between July and August ( Figure 2d; p = 0.005) and a 48-51% increase in 2019 from April to September, only for riparian and grain plots (Figure 2f; riparian p < 0.001, grain p < 0.001, willow p = 0.72). A similar pattern was observed in 2018, a corn year. Willow plots had numerically higher invertebrate diversity and richness in early April compared to the grain and riparian plots. However, over the season, grain and riparian plots had a general increase in both diversity and richness, resulting in significantly lower diversity and richness in willow plots by September in some cases (Figure 2b,e; grain vs. willow diversity p = 0.03; riparian vs. willow diversity p = 0.47; grain vs. willow richness p = 0.001; riparian vs. willow richness p = 0.1). had numerically higher invertebrate diversity and richness in early April compared to the grain and riparian plots. However, over the season, grain and riparian plots had a general increase in both diversity and richness, resulting in significantly lower diversity and richness in willow plots by September in some cases (Figures 2b and 2e; grain vs. willow diversity p = 0.03; riparian vs. willow diversity p = 0.47; grain vs. willow richness p = 0.001; riparian vs. willow richness p = 0.1). The sampling date and seasonality were also found to influence invertebrate community composition more than the landcover type. The PCA biplots ( Figure 3) across all three years show some grouping by sampling date (Figures 3a-3c) and less grouping by landcover type (Figures 3d-3f). With only two sampling dates in 2017, invertebrate community composition (based on species abundance and presence) showed distinct groupings of samples by sampling date, which were associated with 8 invertebrate species in July and 20 species in August. The RF model agreed with a 2% out-of-bag (OOB) error when the sampling date was used as the response variable. With increased sampling in 2018 and 2019, the RF models found that neither sampling date nor landcover type could explain invertebrate community composition alone, with OOB errors ranging from 34% to 60%. The sampling date and seasonality were also found to influence invertebrate community composition more than the landcover type. The PCA biplots ( Figure 3) across all three years show some grouping by sampling date (Figure 3a-c) and less grouping by landcover type (Figure 3d-f). With only two sampling dates in 2017, invertebrate community composition (based on species abundance and presence) showed distinct groupings of samples by sampling date, which were associated with 8 invertebrate species in July and 20 species in August. The RF model agreed with a 2% out-of-bag (OOB) error when the sampling date was used as the response variable. With increased sampling in 2018 and 2019, the RF models found that neither sampling date nor landcover type could explain invertebrate community composition alone, with OOB errors ranging from 34% to 60%.
In both years, the RF model had greater accuracy in predicting invertebrate communities present in early April (OOB errors 10-28%) or in September of 2018 (13% error) than many of the other sampling dates (errors > 46% for both years). The k-means cluster analysis showed similar results, with seasonality (early or late season) in combination with some landcover effect, influencing cluster combinations for all three years (Table 1).
Additionally, the first two PC axes created from the cluster analysis for all years explained a large portion of the variation in invertebrate community composition between sampling locations (35-50% of variation explained).
The three landcover types were found to support a number of unique species (species only observed in a single landcover during the study period), including 21 species (14% of the species) in grain plots, 15 species (10%) in willow plots, and 11 species (7%) in riparian plots. Many of these species were rare and observed only once during sampling. The unique species that we observed at least more than once in each landcover are listed in Table 2. In both years, the RF model had greater accuracy in predicting invertebrate communities present in early April (OOB errors 10-28%) or in September of 2018 (13% error) than many of the other sampling dates (errors > 46% for both years). The k-means cluster analysis showed similar results, with seasonality (early or late season) in combination with some landcover effect, influencing cluster combinations for all three years (Table 1). Table 1. Results of K-means cluster analysis for the invertebrate community broken down by sampling year, number of clusters identified, number of samples from each landcover type, number of invertebrate species associated with each cluster, and a description of the make-up of each cluster by landcover, sampling date, and most abundant invertebrate species.   In contrast, 42% of the 155 species observed were found commonly across all three landcovers (Supplementary Materials Table S1). An additional 16% of species were common between willow and grain plots, compared to only an additional 8% that were common between grain and riparian plots or the 4% that were common between riparian and willow plots. There were 14 different invertebrate orders observed across the study period ( Figure 4), with new orders observed each year.

Year
For example, the orders Entomobryomorpha (springtails: associated with willow and riparian plots) and Ixodida (ticks: associated with willow and grain plots) were new orders observed in 2019. A potential crop effect was also observed, where orders Orthoptera (grasshoppers/crickets) and Odonata (dragonflies/damselflies) were only observed during soybean production years (2017 and 2019). There was also an increase in the abundance of pest species, including Thysanoptera (Thrips palmi), which was associated with all three landcovers in 2019, but the largest abundances were seen in soybean plots. The invertebrate community observed across landcover types and sampling years also varied widely in the functional roles of species present ( Figure 5). In 2017 and 2018, willow plots supported the highest number of species across functional groups in four of the six sampling periods. Willow plots supported a larger number of pollinators in the spring of 2018, although the number of pollinators that was caught at that time was relatively low. In 2019, the number of observed pollinators across landcovers increased, with soybean plots found to support a greater number. In fact, across the entire 2019 season, soybean plots consistently supported the greatest number of invertebrate species, with 50 different species observed across soybean plots in September alone. The majority of these species were classified as either pest, pollinator, or predator. The classification of species by functional group is listed in Supplementary Materials Table S1. In contrast, 42% of the 155 species observed were found commonly across all three landcovers (Supplementary Materials Table S1). An additional 16% of species were common between willow and grain plots, compared to only an additional 8% that were common between grain and riparian plots or the 4% that were common between riparian and willow plots. There were 14 different invertebrate orders observed across the study period ( Figure 4), with new orders observed each year.   For example, the orders Entomobryomorpha (springtails: associated with willow and riparian plots) and Ixodida (ticks: associated with willow and grain plots) were new orders observed in 2019. A potential crop effect was also observed, where orders Orthoptera (grasshoppers/crickets) and Odonata (dragonflies/damselflies) were only observed during soybean production years (2017 and 2019). There was also an increase in the abundance of pest species, including Thysanoptera (Thrips palmi), which was associated with all three landcovers in 2019, but the largest abundances were seen in soybean plots. The invertebrate community observed across landcover types and sampling years also varied widely in the functional roles of species present ( Figure 5). In 2017 and 2018, willow plots supported the highest number of species across functional groups in four of the six sampling periods. Willow plots supported a larger number of pollinators in the spring of 2018, although the number of pollinators that was caught at that time was relatively low. In 2019, the number of observed pollinators across landcovers increased, with soybean plots found to support a greater number. In fact, across the entire 2019 season, soybean plots consistently supported the greatest number of invertebrate species, with 50 different species observed across soybean plots in September alone. The majority of these species were classified as either pest, pollinator, or predator. The classification of species by functional group is listed in Supplementary Materials Table S1.

Understory Plant Community
In contrast to the invertebrate results, landcover type had a large effect on the understory plant diversity in both sampling years, 2018 and 2019 (p < 0.001). In 2018, willow plots had a significantly higher chance of having plant diversity greater than zero compared to grain plots, with the chance of observing quadrats without any understory spe-

Understory Plant Community
In contrast to the invertebrate results, landcover type had a large effect on the understory plant diversity in both sampling years, 2018 and 2019 (p < 0.001). In 2018, willow plots had a significantly higher chance of having plant diversity greater than zero compared to grain plots, with the chance of observing quadrats without any understory species decreasing through the season (Figure 6a; p < 0.001). Riparian plots, in contrast, always had plant species present. When the linear model was rerun without absences (0 s) included, willow and grain plots had similar plant diversity (p = 0.25), while they both had significantly lower diversity than riparian plots (Figure 6b; p < 0.001). Similar results were observed in species richness, with a significantly greater number of species found in the riparian plots compared to the willow and corn plots in 2018 (Figure 6d,e; p < 0.001).  The increased sampling in 2019 resulted in fewer observed locations without understory plant coverage (absences). In contrast to the 2018 results, riparian plots only had higher plant species diversity than willow and grain plots until September of 2019 (p < 0.003 for all comparisons), when willow and riparian plot diversity were found to be statistically similar (Figure 6c; p = 0.97). Willow plot diversity significantly increased by 49% across the season (p < 0.0001), with a similar trend observed in species richness through a 42% increase from April to September (Figure 6f; p < 0.002 for all comparisons). Soybean plots, for both plant diversity and richness, had a general decline over time in 2019, with 53% fewer species at the end of the season compared to the spring (p < 0.002 between April and September) and a corresponding 49% reduction in plant diversity (p < 0.04).
Landcover type also played a large role in the observed understory plant community composition. The PCA biplots (Figure 7) show that the first two PCA axes distinctly separate the plant communities by landcover type in both years (Figures 7a and 7b). The re- The increased sampling in 2019 resulted in fewer observed locations without understory plant coverage (absences). In contrast to the 2018 results, riparian plots only had higher plant species diversity than willow and grain plots until September of 2019 (p < 0.003 for all comparisons), when willow and riparian plot diversity were found to be statistically similar (Figure 6c; p = 0.97). Willow plot diversity significantly increased by 49% across the season (p < 0.0001), with a similar trend observed in species richness through a 42% increase from April to September (Figure 6f; p < 0.002 for all comparisons). Soybean plots, for both plant diversity and richness, had a general decline over time in 2019, with 53% fewer species at the end of the season compared to the spring (p < 0.002 between April and September) and a corresponding 49% reduction in plant diversity (p < 0.04).
Landcover type also played a large role in the observed understory plant community composition. The PCA biplots (Figure 7) show that the first two PCA axes distinctly separate the plant communities by landcover type in both years (Figure 7a,b). The results of the RF model also confirm these results, with OOB errors of less than 3% in predicting sample origin by landcover type based on the plant community present. In both years, red-seeded dandelion (Taraxacum erythrospermum), which was highly abundant in willow plots, was the largest predictor in the models for sample discrimination. Common blue violet (Viola sororia) was another top predictor for both years, along with perennial ryegrass (Lolium perenne) and poison hemlock (Conium maculatum) in 2018 and Virginia ryegrass (Elymus virginicus) in 2019. The results of the PCA analysis ( Figures  7c and 7d) indicated that the sampling date and seasonality also play roles in understory plant community composition, which is supported by the k-means cluster analysis (Table  3). Most clusters were comprised of samples from the same landcover type, but landcovers were found to be split by species composition during different parts of the season, as well as landscape position (marginal versus non-marginal soils). Grain plots, in particular, were more clearly divided by seasonality, with different communities present in the spring versus the fall.
Of the 77 plant species observed in the 2 sampling years (Supplementary Materials  Table S2), the willow plots had the highest number of total species observed (56), followed by the riparian (52) and grain plots (24 total species). The riparian and willow plots had the greatest number of plant species in common with one another (25), while another 13 plant species were found to be common across all three landcovers. As observed in the invertebrate community, each landcover also supported unique species (willow-16 ; riparian-12; grain-7). Unique species that were observed more commonly under each landcover type are shown in Table 2. Between the two sampling years, the majority of understory plants were native species across all three landcovers (grain: 50%; willow: 58%; riparian: 73%). For riparian and grain plots, native species tended to have greater relative ground coverage as well (a measure of abundance) compared to introduced spe- In both years, red-seeded dandelion (Taraxacum erythrospermum), which was highly abundant in willow plots, was the largest predictor in the models for sample discrimination. Common blue violet (Viola sororia) was another top predictor for both years, along with perennial ryegrass (Lolium perenne) and poison hemlock (Conium maculatum) in 2018 and Virginia ryegrass (Elymus virginicus) in 2019. The results of the PCA analysis (Figure 7c,d) indicated that the sampling date and seasonality also play roles in understory plant community composition, which is supported by the k-means cluster analysis (Table 3). Most clusters were comprised of samples from the same landcover type, but landcovers were found to be split by species composition during different parts of the season, as well as landscape position (marginal versus non-marginal soils). Grain plots, in particular, were more clearly divided by seasonality, with different communities present in the spring versus the fall.
Of the 77 plant species observed in the 2 sampling years (Supplementary Materials  Table S2), the willow plots had the highest number of total species observed (56), followed by the riparian (52) and grain plots (24 total species). The riparian and willow plots had the greatest number of plant species in common with one another (25), while another 13 plant species were found to be common across all three landcovers. As observed in the invertebrate community, each landcover also supported unique species (willow-16; riparian-12; grain-7). Unique species that were observed more commonly under each landcover type are shown in Table 2. Between the two sampling years, the majority of understory plants were native species across all three landcovers (grain: 50%; willow: 58%; riparian: 73%). For riparian and grain plots, native species tended to have greater relative ground coverage as well (a measure of abundance) compared to introduced species (Figure 8). For willow plots, however, introduced species were found to have greater ground coverage than native species, which was primarily due to the large presence of red-seeded dandelion (Taraxacum erythrospermum) in the northern willow plots.
In the willow and riparian plots, the majority of vegetation species were perennials, while annuals were more common in the grain plots. Additionally, the majority of species observed across all landcover types tended to be forbs (broad-leaf). Forb species also had greater ground coverage in both years in willow plots, while it varied for grain (vine or forbs) and riparian (tree or graminoid) plots depending on the year.

Vegetation Community Influence on Invertebrate Community
A symmetric co-correspondence analysis (Co-CA) was used to assess the covariation/commonality between the vegetation and invertebrate community matrices. The Co-CA axes from the 2018 data were found to explain 71% of the total variation in the invertebrate matrix and 92% of the total variation in the vegetation matrix. The trend was the same for 2019, with 79% of the variation explained in the invertebrate community and 97% in the vegetation community. The biplots in Figure 9 show the commonality between the two matrices. Table 3. Results of K-means cluster analysis for understory vegetation community broken down by sampling year, number of clusters identified, number of samples from each landcover type, number of plant species associated with each cluster, and a description of the make-up of each cluster by landcover, sampling date, and species with the higher relative percent cover.

Vegetation Community Influence on Invertebrate Community
A symmetric co-correspondence analysis (Co-CA) was used to assess the covariation/commonality between the vegetation and invertebrate community matrices. The Co-CA axes from the 2018 data were found to explain 71% of the total variation in the invertebrate matrix and 92% of the total variation in the vegetation matrix. The trend was the same for 2019, with 79% of the variation explained in the invertebrate community and 97% in the vegetation community. The biplots in Figure 9 show the commonality between the two matrices. Both the influence of the sampling period and the landcover type could be observed in the grouping of invertebrate and vegetation communities. As observed in PCAs, willow and grain vegetation communities tended to cluster separately. Grain vegetation communities were also observed to differ by sampling period with two groups (April and June-September). Willow vegetation communities tended to be clustered together regardless of sampling period across the two years. For the invertebrate communities, April communities tended to group together in both years. In 2018, the Co-CA2 axis was observed to separate the April grain invertebrate communities (more positive Co-CA2 scores) and the April willow and riparian invertebrate communities (negative Co-CA2 scores).

Vegetation Community Influence on Invertebrate Community
A symmetric co-correspondence analysis (Co-CA) was used to assess the covariation/commonality between the vegetation and invertebrate community matrices. The Co-CA axes from the 2018 data were found to explain 71% of the total variation in the invertebrate matrix and 92% of the total variation in the vegetation matrix. The trend was the same for 2019, with 79% of the variation explained in the invertebrate community and 97% in the vegetation community. The biplots in Figure 9 show the commonality between the two matrices. Both the influence of the sampling period and the landcover type could be observed in the grouping of invertebrate and vegetation communities. As observed in PCAs, willow and grain vegetation communities tended to cluster separately. Grain vegetation communities were also observed to differ by sampling period with two groups (April and June-September). Willow vegetation communities tended to be clustered together regardless of sampling period across the two years. For the invertebrate communities, April communities tended to group together in both years. In 2018, the Co-CA2 axis was observed to separate the April grain invertebrate communities (more positive Co-CA2 scores) and the April willow and riparian invertebrate communities (negative Co-CA2 scores). Both the influence of the sampling period and the landcover type could be observed in the grouping of invertebrate and vegetation communities. As observed in PCAs, willow and grain vegetation communities tended to cluster separately. Grain vegetation communities were also observed to differ by sampling period with two groups (April and June-September). Willow vegetation communities tended to be clustered together regardless of sampling period across the two years. For the invertebrate communities, April communities tended to group together in both years. In 2018, the Co-CA2 axis was observed to separate the April grain invertebrate communities (more positive Co-CA2 scores) and the April willow and riparian invertebrate communities (negative Co-CA2 scores).
The influence of the vegetation community composition on the invertebrate community composition was evaluated with a predictive Co-CA, which used the vegetation community matrix as a predictor. None of the Co-CA axes in either year were found to be significant based on the cross-validation method. Overall, the results from both years found that the understory vegetation community was not a good predictor for the invertebrate community composition.

Discussion
A total of 155 invertebrate species and 77 understory plant species were observed during this study across all landcover types and sampling years (Supplementary Materials  Tables S1 and S2). The influence of landcover type on invertebrate and plant diversity and richness varied. Landcover type was found to influence species diversity and richness more for understory plants than invertebrates, but seasonality played a larger role for invertebrates. For vegetation, as expected, the riparian plots had significantly higher species diversity and richness than the willow and grain plots for the majority of two sampling seasons. However, this was not always true for invertebrates. When statistical differences were observed in invertebrate diversity or richness among landcover types, generally grain and riparian plots were similar to each other and had higher diversity and richness than observed in willow plots. This differed from previous studies that found higher arthropod diversity and activity in willow than in other arable crops, such as corn and break crops [12,43]. Verheyen et al. [12] found several invertebrate orders (Hymenoptera (bees, ants, wasps); Coleoptera (beetles); Araneae (spiders); and Diptera (flies)) to have higher abundances in SRWC plantations and a greater number of invertebrates among functional groups (omnivores, detritivores, mycophages, phytophages, and parasitoids) than what was found in corn plots. In this study, however, the majority of the aforementioned orders and functional group abundance (pests, predators, and beneficial) were found to be rather similar between willow and grain plots across the three years. The exception was Araneae, which was also found to be in higher abundance in willow and riparian plots compared to grain plots. However, other sampling techniques, such as pitfall traps, may be better suited to monitoring some of these invertebrate orders.
Visual differences in invertebrate diversity and richness were also seen between the two grain types (corn and soybean). The soybean years showed higher trends in invertebrate diversity and richness than the corn year. Differences in species' order abundance between the years also suggest a crop effect. This differs from the results in Wheelock [16], which observed a greater number of insects in cornfields versus soybean fields in a two-year study, as well as similar communities between the two crops. Wheelock [16] also noted, however, that trap height played a large role in the species abundance and composition observed in the traps. In the case of pollinators, the authors observed a lower abundance of pollinators at ground level and corn ear height compared to traps set at tassel height. In this study, all traps were set at a single height for consistency between crops and sampling years; however, as a result, trap height varied relative to plant height and flowering height. In the case for both willow and corn, trap height was lower than catkin or tassel height, whereas trap height was around flowering height for soybean. Observations during the study commonly found long-horned bees (Melissodes spp.) foraging on pollen from corn tassels, and hoverflies (Syrphidae family) foraging on willow catkins. Even in the spring, in grain plots, when flowering chickweed (Stellaria media) and hembit (Lamium amplexicaule) were abundant, it was observed that fewer pollinators were trapped during this period than were observed foraging near the ground. Trap height may also affect the presence of other functional groups, as in the case for willow, even leaf development starts from the top of the branches, in the upper canopy, with lower canopy leaves developing later in the season. Designing traps that include bee bowls placed at multiple heights as discussed in Wheelock [16], or using other trapping techniques such as sweep netting, may be able increase the invertebrate community sampled.
Trap height may have also impacted conclusions on whether willow has served as an important food source in the early spring for pollinators, which was reported by other studies [8,11,13,18]. Numerically, willow plots had higher diversity and richness in early April of 2018 compared to riparian and grain plots, as well as a greater number of observed pollinators. However, in both cases, this trend was not observed in 2019, a soybean year. This was initially surprising, as the visual number of catkins observed was higher in the spring of 2019 than in 2018, a factor of the short harvest cycle. Willows are mainly insectpollinated (entomophilous), and some bee species in North America use willow as their main pollen source, including many Andrena species, an important pollinator for berry and fruit crops [13,44]. Only one Andrena species was observed during the study, Andrena nasonii, and only during the soybean years. However, Ostaff et al. [13] found that Andrena species and hoverflies prefer male flowers (which provide pollen and nectar) over female flowers (nectar only). Pollen is a nutritious protein source for growth and reproduction, whereas nectar is an easily digestible energy source [44]. The willow clone used in this study was female. Reddersen [44] also noted that flowering abundance can differ between willow species and clones, highlighting the importance of species/clone selection for pollinators. Volk et al. [8] suggested diversifying species, hybrids, or clones to increase heterogeneity; however, proper selection is equally important to ensure compatibility between willow types so that one does not outcompete the others. The importance of pollinator support in the Midwest U.S. is highlighted by the potential benefits pollinators, particularly native bee communities, can have on soybean seed yields. Although yield benefits are not always seen, Cunningham-Minnick et al. [45] argues that production areas would not receive the benefit of pollinator services without the presence of nearby habitat providing food and nesting areas. Shrub willow has the potential to provide nesting habitat as demonstrated in a modeling study by Graham et al. [18]; however, the authors did find less nesting habitat potential under willow than for other landcover types, including switchgrass and native prairie, arguing the need for additional field-based studies. The presence of understory species that flower at different times during the season may provide additional resources for pollinators and alter nesting potential [18]. Furthermore, an in-field buffer design may also increase the crop-edge habitat that native pollinators may forage along, increasing the potential benefits other field crops may receive [45].
The presence of other understory plant species across landcover type may also be important in supporting invertebrates and pollinators throughout the season. Although the Co-CA did not find the understory vegetation community to be a good predictor for invertebrate community composition, there was a high amount of covariation between the communities. Furthermore, the presence of species such as Halictinae (Agapostemon sweat bee), a pollinator that is primarily active later in the season, in willow plots suggests the understory plant community or general conditions in the willow plots may help support species throughout the season. The presence of winter weeds (chickweed and henbit) in the spring in grain plots were also found to support a number of different pollinators. Furthermore, the use of willow edges (headlands) can be another opportunity to increase plant diversity of flowering species. In this study, several plant species, including redseeded dandelion (Taraxacum erythrospermum), common blue violet (Viola sororia), goldenrod species (Solidago spp.), and daisy fleabane (Erigeron strigosus), were highly attractive to a range of foraging invertebrate orders, with greater abundances of these plant species observed along the edges of the willow plots. Willow catkins in the spring were also observed to be more abundant along the edges of the plots or in less dense areas of the plots, in addition to the upper canopy. These areas were not well-represented during sampling, as trap placements in willow plots were well within the borders of the plots.
This highlights, however, the potential benefits of bioenergy landscape designs, such as buffers, which have a higher proportion of edge area than a dedicated field. Previous research has shown that headlands and field margins are important for predatory invertebrate species, such as beetles and spiders, as they provide refuge during the winter as well as corridors for movement during the year [46]. The importance of predator presence for biological control is highlighted by the general increase in pest species observed across the seasons in this study. Thrips palmi, an agricultural pest that feeds on the leaves, stems, and flowers/fruits of a wide range of hosts [47], was commonly found across all landcover types, with the greatest abundance in soybean plots in 2019. However, natural predators of Thrips were also present across the three landcovers, including the orders Neuroptera (lacewing flies) and Hemiptera (true bugs), such as the minute pirate bugs (Miridae) [48][49][50]. Minute pirate bugs find willow, among other plant species, attractive habitat [49]. There were also a wide range of predators found across all three landcovers. Araneae (spiders) and adults of Coccinellidae (ladybugs) can help regulate aphids and other arthropod populations, and Carabids, which were commonly found, help regulate pests along with weed species [51,52]. Furthermore, previous research has found edges of SRWC to be favored by other trophic levels, such as birds and mammals, resulting in suggestions that SRWC production designs that are smaller and have higher edge-to-area ratios are preferred to increase wildlife and invertebrate benefits [46]. With the increased introduction of bioenergy crops and the scale-up of production, Landis and Werling [53] do argue that positive and negative changes in the invertebrate communities could arise. They highlighted an example when corn production increased with the demand for ethanol, resulting in soybean yield loss and added pest management costs (USD 58 million per year) due to the reduction in populations of a natural aphid predator. Pre-and post-bioenergy crop introduction studies may be helpful in addressing this uncertainty and determine if the introduction of undesirable species, that can harm surrounding crops, occurs. Additionally, shrub willow, similar to most crops, can be susceptible to invertebrate pests that can impact biomass productivity. Many management approaches have already been recommended to address this risk including cultivar selection and planting smaller areas instead of plantations [53]. Therefore, integrated bioenergy buffers, such as the design used in this study, may have reduced pest risk potential.
The benefit of willow production in this study was also highlighted by the presence of unique invertebrate and understory plant species. Willow plots supported nine unique plant species and three unique invertebrate species (counting species that were observed more than once). When rare invertebrates are included, the majority of unique invertebrates found within the willow plots were classified as predators. For plants, the majority of the unique species were perennial, native species. The majority of plant species observed in willow plots were native species (58% between the two years). However, the introduced species tended to have greater ground coverage. This was primarily due to the large presence of red-seeded dandelion in the northern willow plots. As the cluster and random forest analyses highlighted, the presence and abundance of red-seeded dandelion was a large contributing factor to the distinction between understory plant communities of willow plots and either riparian or grain plots in this study.
In grain plots, annual species were found to be more common and more abundant than perennial species, a contrast to the community composition of willow and riparian plant communities. The greater presence of annual species, as well as the lower diversity and richness of plants within grain plots, is likely due to herbicide and tillage treatments, which do not occur in willow or riparian areas. The greater presence of perennial species in willow plots does, however, support previous work that found an increase in perennial species post-SRWC establishment [46]. Willow plots also had comparable species diversity and richness to riparian plots by the end of the 2019, although the density or ground coverage of each species was generally far less than that observed of species found in the riparian plots. Riparian and willow plots also had a larger number of plant species in common than those found in grain plots. This differs from what was observed in the invertebrate community, which found willow plots to have more species in common with grain plots than with riparian plots.

Conclusions
This 3-year study aimed to evaluate an additional non-target benefit (supporting biodiversity) of strategically placed short-rotation shrub willow bioenergy buffers within a Midwest U.S. grain field, buffers that were initially designed for addressing nutrient loss reduction. The study took place in Fairbury, IL during the second willow harvest cycle (2017-2019). The canopy invertebrate community was monitored using bee bowls, and the understory (non-target crop) plant community was assessed using point-sampling ground coverage measurements. The results of this study showed that the agricultural landscape system, a combination of grain crop (corn/soybean; managed with herbicide and/or tillage), shrub willow (minimally managed for non-target species), and riparian-edge habitat, supported a large number of invertebrate (155) and plant (77) species. Landcover type had a greater impact on understory plant community composition, diversity, and richness, with distinct communities observed across all three landcovers. Riparian plots had greater plant species diversity and richness than willow or grain plots for most of 2018 and 2019. However, by the end of 2019, the willow plots had significantly greater plant diversity and richness than the grain plots and were comparable to the riparian plots. In contrast, the willow plots generally had lower invertebrate diversity and richness than the grain or riparian areas. Trap height could have affected these results. Seasonality had a larger influence on diversity and invertebrate community composition than landcover type, with the majority of species found commonly among the grain, willow, and riparian plots. Unique plant and invertebrate species, however, were observed across all three landcovers, highlighting the importance of each habitat for supporting biodiversity. Although, the buffer placement design in this study targeted nutrient loss as the primary ES, the presented results and other studies found in literature highlight the potential for bioenergy landscape designs, such as buffers, for supporting biodiversity and thereby providing an additional ecosystem service by an integrated perennial bioenergy crop.