Changes in conservation value from grasslands to savannas to forests: How a temperate canopy cover gradient affects butterfly community composition

Temperate savannas and grasslands are globally threatened. In the Midwest United States of America (USA), for example, oak savannas persist today at a small percentage of recent historic coverage. Therefore, restoration of habitats of low and intermediate canopy cover is a landscape conservation priority that often emphasizes returning tree density to a savanna-like target value. Understanding how animal species react to such changes in vegetation structure is important for assessing the value of these restoration plans. We examined how butterfly community attributes in northwest Indiana USA, including community composition, richness, and abundance responded to a grassland-to-forest gradient of canopy cover. Butterfly community composition under intermediate canopy cover differed significantly from community composition in the most open or closed-canopy habitats. Composition of the plant community in flower was a significant predictor of three assessed attributes of the butterfly community—composition, richness, and abundance. Phenology, expressed as day-of-the-year, was also a strong predictor of these butterfly community attributes. Few butterfly species were habitat specialists as adults although canopy cover was a more important predictor of adult community composition than of richness or abundance of butterflies. Therefore, adult butterfly community differences along the canopy cover gradient were less about butterfly communities filled with habitat specialists for different canopy-defined habitats and more about gradual changes in community composition along this gradient. Overall, butterfly community richness was predicted to peak at about 34% canopy cover, butterfly abundance at about 53% canopy cover, community conservation value at about 59% canopy cover, and a combination of desirable conservation attributes–high diversity, high abundance, and high conservation value–was predicted to reach a peak of co-occurrence at about 67% canopy cover suggesting that habitats of intermediate canopy cover might be particularly effective for butterfly conservation in this region.

Introduction Globally, temperate savannas, woodlands, and grasslands are biomes of special conservation concern [1]. Human activities have influenced the temperate forests, savannas, and grasslands of eastern North America for thousands of years [2,3]. Nonetheless, since settlement by European-Americans, the eastern North American landscape has changed greatly in part due to conversion to agriculture and other development and due to changes in disturbance regimes, especially fire. Additionally, this landscape has experienced relatively low rates of conservation protection [1,4]. Savanna was present on more than ten million hectares (ha) in the Midwest USA and Canada in the early nineteenth century but the rate of conversion to human use and the frequency of succession to more forested habitats, due to decreased fire frequency, have been high even compared to global trends [5][6][7]. Estimates of the decrease in acreage of oak dominated woodlands and savannas in this region range from 83% [8] to more than 99% [7] over the past two hundred years. Recognition of this decline has stimulated planning for large scale restoration of these habitats in the Midwest United States [8,9].
Prairie, savanna, and forest habitats co-occur in a Midwest USA habitat transition zone with the distribution among these habitats across landscapes likely related to climatic factors, such as precipitation levels, to soil types, to landscape physiognomy, and to fire frequency [10]. In deciding habitats to prioritize during restoration along this canopy cover gradient, an understanding of how plant and animal communities are likely to change with changes in canopy cover is valuable. Previously, to provide guidance to managers as they formulate goals for restoration within this habitat mosaic landscape, we examined the effects of canopy cover on amphibians, bees, birds, and reptiles along this gradient, a group of animal taxa that differ in their vagility, size, life span, and life history habitat requirements [11][12][13][14]. Here, we add adult butterflies to this group of studied taxa for several reasons including their complex life cycles that might promote dependencies on multiple canopy conditions. Also, at the time of this study, canopy reduction to provide suitable habitat for the federally endangered Karner blue butterfly (Lycaeides melissa samuelis) was a primary goal for managing some of the sites we studied but the effects of this management goal on other species of butterflies or other animal or plant species was not well documented [15][16][17].
Concern about how changes in canopy cover and heterogeneity, due to changes in fire or cutting regimes, might affect understory communities are geographically widespread [18]. For example, canopy management via coppicing, a woodland management technique historically used in Europe but less consistently used recently, has resulted in significant changes in canopy structure in many areas of Europe [19,20]. These changes in canopy cover have affected the abundance and distribution of a variety of plant and animal species, including arthropods such as spiders [21], beetles [22,23], moths [24], and butterflies [25,26]. Frequently, loss of canopy heterogeneity is associated with decline in light penetration and microclimatic heterogeneity in the understory and decline in species diversity [21]. Here we expand our investigation of how canopy cover manipulation might affect biotic diversity in the Midwest USA prairie-forest transition zone [27] by examining how canopy cover affects adult butterfly community composition, richness, and abundance. We present predictions for the response of these butterfly community attributes to canopy cover. Finally, we predict whether canopy cover might also affect the conservation value, or ability to retain threatened species, at these sites.

Methods
We examined relationships between butterfly distribution and habitat structure along a canopy cover gradient, using an historic data set of butterfly distribution collected in 1998-1999, at 25 sites within three locations in northwest Indiana, USA. Data used for this study are publicly available at https://doi.org/10.5066/P9TYX2AG.

Study sites
The 25 study sites were located at Indiana Dunes National Lakeshore (now National Park) (n = 17 sites) (41˚38' N, 87˚09' W) (6000 ha), Tefft Savanna Nature Preserve and Jasper-Pulaski Fish and Wildlife Area (n = 7 sites) (41˚10' N, 86˚58' W) (3250 ha), and Hoosier Prairie Nature Preserve (n = 1 site) (41˚31' N, 87˚27' W) (225 ha). Sites, and the habitat assessment methodology, are described in detail elsewhere but are reviewed here [11][12][13]28]. At each site, butterflies were surveyed along either one 500-m transect, or two 250-m transects separated by 30 m to limit double counting, depending on site dimensions. Mean site area was 11.8 ha ± 2.4 (standard error, SE). Sites were situated from 0.8-80 km inland from the southern shore of Lake Michigan, averaged 1.8 km ± 0.7 (SE) between nearest neighbors, were mainly on sandy soils, and represented five replicates of each of five habitat types that may arise from historic changes to savannas. These habitat types included: (1) open habitats (< 20% canopy cover) including prairie remnants, other grasslands, and sites with low canopy cover due to past soil disturbances, such as farming and sand mining [29], (2) savannas with 20-50% oak-dominated canopy cover and little woody understory vegetation, (3) woodlands with 50-90% canopy cover and little woody understory vegetation, (4) scrub habitats with high density (> 1,000 woody stems) of black oak (Quercus velutina) sprouts (2.5-10 cm diameter breast height (dbh) ha -1 and < 5 m high) that often arise following intense fires, and (5) forests with high canopy cover (>90%) from larger trees and with several woody vegetation layers containing multiple tree species with more than 300 woody stems > 10 cm dbh ha -1 . Dominant tree species included black oak and white oak (Q. alba) and sassafras (Sassafras albidum).

Butterfly abundance determination
We surveyed butterflies 21 times (n = 13 in 1998 and 8 in 1999) from May 4 -October 13 in 1998 and from April 26 -October 1 in 1999. At least one genus (Phyciodes) contained more than one species that could not be distinguished in the field. In summing species numbers, we counted this genus as one species.
To reduce bias related to differences in detectability among butterfly species, we used distance sampling to convert raw butterfly counts to estimates of density [30,31]. To implement distance sampling, we estimated the perpendicular distance from each observed butterfly to the transect line. If observations of a species were insufficient to estimate the detection function used in the density calculation, we pooled observations over species similar to this rarer species [31]. The distance sampling software [30] calculated effective strip width (ESW), the distance at which we expect the number of a butterfly species detected outside the ESW should equal the number of animals missed inside the ESW. Approximately 7.8% of butterflies observed could not be identified to species in the field but only to a higher taxonomic level, such as family. We apportioned those butterflies to the species level, based on the butterflies identified to species at that site daily, if possible. Density was expressed as the sum of individual species' densities per hectare, either on a per-survey basis or the average per-survey density across all surveys. additional information on a few of the predictors that we do not describe in detail elsewhere [11][12][13]28] or for which the brief description in Table 1 might be insufficient. SEVM: we used an eigenvector based spatial filtering method, SEVM (spatial eigenvector mapping), to account for spatial trends in dependent variables in the analytical models [34,35]. The SEVM filter helps account for the effect of spatial autocorrelation from the relationship between the environmental predictors and butterfly responses.
PlantComposition: in both late summer 1999 and spring/early summer 2000, we sampled vegetation twice in twenty randomly placed 1 x 2 m plots, each containing five 2:1 nested subplots. The smallest subplot in which a plant species was observed was recorded for each plot and a frequency score was assigned, with 5 representing plants found in the smallest subplot and 1 in the largest subplot and 0 when the plant species was not found in the plot [36]. We averaged the late summer and the spring/early summer scores per species across a site separately and then used the greater of the summer or spring scores to represent the maximal abundance of that plant species at that site. These maximal species abundances were then ordinated in a one-dimensional principal curve analysis [37,38]. The ordination score characterized the plant community occurring at a site and was used as a possible predictor of the butterfly community.
FlowerComposition: the total number of stems (ha -1 ) with flowers within 1 m of five counting points along the 500 m transect was counted and a one-dimensional principal curve VegetationShort Vertical percent cover of vegetation < 0.3 m tall [11][12][13]28] VegetationTall Vertical percent cover of vegetation 0.3-1.0 m tall [11][12][13]28] Litter % litter cover [11][12][13]28] CanopyCover % canopy cover [11] SEVM Spatial eigenvector mapping [34,35] Fire2 total area burned over the two years preceding surveys within 200 m of the transect, divided by the total area within 200 m of the transect [28] Fire15 total area burned over the fifteen years preceding surveys within 200 m of the transect, divided by the total area within 200 m of the transect [28] PlantComposition Principal curve ordination score describing the yearly plant community at a site [36][37][38] FlowerComposition Principal curve ordination score describing the plants in flower during a survey [37,38] FlowerSpecies Number of species in flower, estimated by Abundance Coverage-based Estimator [39,40]  ordination score calculated to describe the flowering plants available during each cycle of observation.
FlowerSpecies: for the flowering plants we counted when surveying a transect, we used the daily stem counts during each of the twenty surveys to calculate an Abundance Coveragebased Estimator (ACE) estimate of richness of plants in flower across our surveys at a site [39,40].
FlowerStems: was the number of stems of flowering plants counted during a survey as described for FlowerComposition. Total area of five points surveyed equals 4 m 2 .
Vegetation and site characteristics, including VegetationShort, VegetationTall, Litter, Cano-pyCover across each site were measured in six 0.05 ha (25.2 m diameter) circular plots and averaged across the site by inverse distance weighting from transect center [41].
Three responses to these environmental variables were examined: (1) Butterfly Species Richness-the number of butterfly species observed along a transect survey, (2) Butterfly Composition-the one dimensional principal curve score [38] of the butterfly community composition based on densities of the species observed, (3) Butterfly Density-sum of densities of individual butterfly species during the survey of a transect (# ha -1 ).

Statistical analysis
We used PERMANOVA [42,43] to assess the significance of differences in butterfly community composition among the five habitat types and between pairs of the habitat types. Significance of differences in abundance, richness, and conservation value across habitats was assessed by a Kruskal-Wallis nonparametric one-way analysis of variance and significance of differences between pairs of habitats by Dunn's multiple comparisons' test [44,45].
Co-correspondence analysis compares ordinations of two communities sampled at the same set of sites [46,47]. Predictive co-correspondence analysis examined whether composition of the plant community, both the overall plant community determined by sampling plots and the flowering plant community determined by counting plants in flower during the butterfly surveys, significantly predicted composition of the butterfly community.
Based on species accumulation curves [48], we estimated whether the number of species captured at a common number of individual butterfly sightings (species richness), or at a common number of sites (species density), was significantly different between pairs of habitats, based on a z-test [49] with p values adjusted for multiple tests, using the Benjamini-Hochberg procedure for multiple testing in R [50,51]. We also calculated an Abundance Coverage-based Estimator (ACE) to estimate number of butterfly species predicted to occur across surveys at a site [39,40].
We used boosted regression tree (BRT) analysis [52,53] to determine how environmental predictors might be related to the three responses, Butterfly Species (richness), Butterfly Composition (community composition), and Butterfly Density (abundance). The implementation of BRT in the R dismo package [52] includes the gbm.step and gbm.simplify procedures that use cross validation to evaluate the number of regression trees and the learning rate for performing the BRT analysis. gbm.simplify then calculates a simplified model that takes into account the tradeoff between model fit and complexity [53]. The relative influence of each variable on the response was scaled to sum to 100, with higher numbers indicating stronger influence. Partial dependence plots were used to illustrate the effects of individual predictors on fitted values of the responses while averaging out the effects of the other predictors. After testing for best parameter settings for the BRT models, we used the following parameter settings in the gbm.step procedure, tree.complexity = 3, learning.rate = 0.01, bag.fraction = 0.75. We assessed the fit of the cross-validated BRT models using correlations between the predicted and observed responses first on a subset of the data that was used for training and then on the remaining data that were used to validate the model derived from the training data [52,53]. Squaring the correlations produces a measure of how well the model results account for variation in the observed results.
We calculated nonlinear regressions [54] based on cubic polynomial or Gaussian 3-parameter peak models, choosing best fit, to produce curves describing the relationships between canopy cover and butterfly community attributes including composition, richness, abundance, and conservation value.
Conservation value. We calculated a butterfly community conservation value based on species' NatureServe conservation status rank score: (1) critically imperiled, (2) imperiled, (3) vulnerable, (4) apparently secure, (5) secure [55]. Prior to weighting, we reversed the scale, so scores ranged from 1 (secure) to 5 (critically imperiled). Higher values then corresponded to more imperiled species. The conservation value index (CVI) was the product of the reversed NatureServe conservation status rank multiplied by a butterfly's density, summed across butterfly species, and was calculated for each survey day for each site. We limited the calculation to those butterflies with reversed NatureServe ranks greater than 2, including only species considered vulnerable, imperiled, or critically imperiled in the index calculation.
Data collection was carried out under collection permits from the US National Park Service, Indiana Dunes National Lakeshore (now Indiana Dunes National Park) and the Indiana Department of Natural Resources, Division of Nature Preserves and endangered species permits from the US Fish and Wildlife Service. All data were collected on public lands managed by the US National Park Service or Indiana Department of Natural Resources.

Results
Sixty-one butterfly species (counting Phyciodes spp. as one species) were observed among the 10,041 butterflies counted in 525 surveys. Mean densities of individual butterfly species ranged from 0.007 to 14.9 ha -1 when averaged across all 21 surveys at all 25 sites (Table 2). Effective strip width (ESW), ranged from 0.96 to 8.5 m (mean 3.24 m ± 0.25 SE), a nearly nine-fold difference in detectability that was used to adjust the area over which the counts were effectively made.

Relationship of butterfly community to plant community composition
Both composition of the overall plant community and composition of the community of plants in flower at the time of survey were significant predictors of butterfly community composition. We performed two sets of co-correspondence analyses (Fig 1). The first analysis examined whether the plant community at a site significantly predicted the butterfly community. In this analysis, the cross-validated fit (8.1%) was nearly at a maximum with two co-correspondence analysis axes. Positive fits are typically considered significant for co-correspondence analyses in which negative fits are possible. The first two axes explained 35.4% of the variance in butterfly community composition accounted for by the 24 axes calculated. The second analysis examined whether the community of plants in flower at the time of survey significantly predicted the butterfly community. In this analysis, the cross-validated fit was nearly at a maximum with four co-correspondence analysis axes at a leave-one-out cross-validated fit of 7.0%. The first four axes explained 37.0% of the variance in butterfly community composition accounted for by the 24 axes.

Predictors of butterfly community composition
One dimensional principal curve ordination scores were calculated to represent plant and butterfly community composition. The principal curve ordinations accounted for (a) 38.0% of the variation in flowering plant community composition (FlowerComposition, ln (X + 0.5) transformed) (n = 525 surveys), (b) 53.2% of the variation in butterfly community composition (Butterfly Composition, butterfly species densities (square root transformed)) (n = 525 surveys); (c) 58.4% of the variation in the plant community composition (PlantComposition) (n = 25 sites).
Butterfly density peaked at about 53% canopy cover while species richness peaked at lower canopy cover, approximately 34% (Fig 2A and 2B). Per-survey density was highest in woodlands and lowest in forests ( Table 3). The relationship between the butterfly community principal curve ordination score (Butterfly Community) and canopy cover (Fig 2C) suggests a gradual change in butterfly community composition from 0% to about 73% canopy cover at sites and then limited butterfly community compositional change at higher canopy covers. Permutational multivariate analysis of variance (PERMANOVA) indicates that these changes in composition with increased canopy cover caused butterfly community composition to differ significantly between open habitats and woodland, scrub, and forest habitats while butterfly community composition in forest habitats was significantly different than in the other four habitat types (Table 4). Canopy cover affected distribution of butterflies by family. Significantly more observations of butterflies in the family Nymphalidae occurred under canopy of trees or shrubs than was true for family Pieridae (Table 5).
Squared cross-validated correlations, ranged from 0.51-0.63 for the three full BRT models describing how well the twenty-one possible variables predicted the three butterfly community responses: composition, richness, and abundance (density) ( Table 6). Except for butterfly species richness, simplified BRT models contained few predictors and all included day-of-theyear (DayofYear) and flowering plant community composition (FlowerComposition). In addition to those two predictors, the simplified model for butterfly community composition included canopy cover (CanopyCover). Twelve predictors were included in the simplified BRT model for butterfly community richness, FlowerComposition, PlantComposition, DayofYear, FlowerStems, HabitatDiversity, Fire15, VegetationShort, CloudCover, CanopyCover, Agriculture, and SEVM in order of importance. Although both composition of the overall plant community and the community of plants in flower were included as predictors in simplified models of at least one of the attributes, richness, composition, or abundance, of butterfly community structure, in general, the composition of plants in flower was a more important predictor of butterfly community characteristics than was the overall composition of the plant community at a site. Based on species accumulation curves (Fig 3A), butterfly species richness was higher in scrub habitats than in the other habitats and species density was highest in open habitats (Fig 3B). Partial dependence plots illustrate the effect of the five predictors with the highest importance values for each butterfly community response after averaging out effects of the other predictors (Fig 4). We define effect size as the range of change in the response (maximum minus minimum) accounted for by a predictor while averaging out effects of the other predictors. For number of species observed per survey (Richness), effect sizes for the five predictors were between 1.42-2.82 species (range of observed values 0-14 species per survey). For community composition (Composition), measured as principal curve ordination value, predictor effect size was between 0.08-0.44 out of a range of possible ordination scores of 0-1. For density (Density), effect sizes were between 59.3-143.9 ha -1 (range of observed values 0-801).
As of 2017, the 61 butterfly species in Table 2 (Table 3). Nonlinear regression of conservation value versus canopy cover suggests that conservation value peaked near 59% canopy cover (Fig 2D).
We examined the covariation of butterfly richness (ACE), density, and conservation value (CVI) per site (n = 25) in Fig 5. Generally, the three conservation responses were not significantly correlated with each other except for density and conservation value (Spearman rank correlation, r s = 0.90, p < 0.001, n = 25, all other pairwise correlations p > 0.05). In Fig 5 we see that the most desirable conservation outcomes-high richness, high density, and high conservation value-did not consistently co-occur, although forested sites were not among the sites with relatively high values of each outcome. We ranked each of the 25 sites from 1 (least desirable outcome: lowest richness, density, CVI) to 25 (most desirable outcome: highest richness, Table 3

. Mean ± SE of butterfly density (ha -1 ) and conservation value index (CVI) of butterfly communities in northwest Indiana by habitat type (n = 5 sites per habitat type).
Within a column, habitats with the same superscript letter are not significantly different from each other (Kruskal-Wallis one-way analysis of variance Dunn's nonparametric multiple comparisons' test [44,45,50] density, CVI) and averaged the three ranks at each site. The averaged rank is plotted against canopy cover in Fig 6. Averaged rank increased from low to intermediate canopy cover (ca. 0-67%) and then decreased. Overall mean is based on all observations within a family. Species mean is the average percentage of observations (± SE), by species, of the # species within a family for which more than twenty observations were made. 1 Family species means followed by the same letter do not differ significantly from each based on Tukey multiple comparisons test [50].
https://doi.org/10.1371/journal.pone.0234139.t005 Table 6. Relative importances for predictors of butterfly species richness, butterfly composition (principal curve score), and butterfly density as determined by boosted regression tree analysis, using procedure gbm.step [52,53]. Importances with daggers ( †) are associated with predictors that were not included in a simplified model (as determined by gbm.simplify). Ten-fold cross-validated correlation between observed and predicted response are shown at the bottom for full and simplified models as well as squared correlation for full model. Model represents predictions for individual surveys (n = 21 surveys per site) across 25 sites (total n = 525).

Discussion
The region of the Midwest USA studied here is within an ecological transition zone where a western grassland dominated ecoregion meets an eastern hardwood forest ecoregion [27]. In this transitional landscape, a diversity of canopy covers occur and manipulation of tree density, using fire, cutting, and planting, is a common form of management [56]. This canopy diversity is likely important for maintaining resident animal and plant diversity [57]. As is true globally,

PLOS ONE
interactions among climate, tree cover, and fire can affect the stability of habitat composition of a landscape that encompasses the canopy cover continuum from grasslands to savannas to forests [58][59][60]. In this Midwest USA region, change in fire regimes, especially fire suppression, has contributed to considerable change in canopy cover and increased local species turnover and has affected regional biodiversity. Ladwig et al. [57], for example, documented 60 years of regional change in plant community composition in a Midwest prairie-savanna-forest mosaic in Wisconsin, USA and found a near doubling of canopy cover in savannas and replacement of prairie-savanna species by forest and non-native species. However, the net effect has often not been a change in richness, but a change in composition. Maintenance of richness but with significant change in composition has been documented in many biodiversity time series globally [61,62]. Canopy cover can affect habitat quality for adult butterflies, and subsequently affect diversity, by modifying local microclimate and light intensity [16,21,26,63,64]. Therefore, understanding responses of vegetation and animals to variation in tree density and shading are important for improving predictability of restoration outcomes and effects of landscape changes if the goal is to maintain historic species composition across heterogeneous landscapes such as the one studied here [18,65]. Indeed, about 8% of the butterfly species from North America north of Mexico (61 of about 765) were found in the relatively small area surveyed in this study, consistent with the great range of canopy cover (~0-100%) that characterized these sites [66,67].
Composition of the butterfly community changed gradually as canopy cover increased from about 0 to 73% cover, then exhibited less community change across higher canopy covers. The result was a forest butterfly community that was significantly different than in lower canopy cover habitats and a butterfly community in sites with the lowest canopy cover (open habitats) that was different than the communities associated with higher canopy cover (woodlands, scrublands, forests). Thus, there were different butterfly community compositions associated with low, intermediate, and high canopy cover. However, change in butterfly community composition was gradual across the canopy cover gradient and there were few butterfly species that were documented as obligates or specialists for the different habitats, open, savanna, woodland, scrub, and forest, defined by canopy cover. Only four of sixty-one butterfly species were classified as habitat specialists and these were associated with low canopy cover. Therefore, gradual compositional changes rather than abrupt turnover of habitat obligate species characterized the relationship between canopy cover and butterfly community composition.
Overall butterfly abundance (density) peaked at intermediate canopy cover, ca. 57%, and the most species (46%) were at their highest concentration in open canopy habitats while the fewest species (5%) were at their highest concentration in forest habitats with the highest canopy cover. Ten (16%) species were most concentrated in savannas, including the Karner blue butterfly, which was the subject of management emphasis at the time of this study but which has since been extirpated from the study area [17]. Richness peaked at lower canopy cover, ca. 34%. Conservation value peaked at intermediate canopy cover, ca. 58%, with the highest mean values in woodland and scrub habitats. Therefore, overall, the butterfly community exhibited gradual change in composition with increasing canopy cover resulting in low, intermediate, and high canopy cover assemblages with few habitat obligate species, overall abundance peaking at intermediate canopy cover but the most species being most concentrated in low canopy cover habitats. Overall richness, measured by species accumulation curves, was similar across habitat types but was predicted to be highest in early-to mid-successional scrub habitats. Species density, or per-survey number of species, was predicted to be highest in open habitats. These patterns of richness and abundance are similar to patterns observed elsewhere. Butterfly community richness and abundance increased in Central European woodlands [26] and Japanese forests and plantations [68] across early successional stages of forest development spanning the most open to low canopy cover stages. That pattern was true for resident, but not migrant, species in the Central European study. For all species in that study, richness and abundance declined with succession from mid-successional to forested stages, much as we documented an increase in richness and abundance from low to mid-level canopy cover and then decreases at higher canopy cover in the current study.
Patterns of variation in butterfly community composition, richness, and abundance were underlain by phenological variation. Day-of-year and composition of the community of plants in flower at the time of survey were among the most important predictors of each of those butterfly community attributes in simplified predictive models. Richness peaked in mid-season and plateaued between 500-1000 flowering stems per 4 m 2 , the total area sampled for flower abundance at each transect. Beyond patterns of seasonal phenological variation, however, butterfly community composition, richness, and abundance were best predicted by different suites of environmental attributes. In particular, the simplified model describing butterfly community richness had twelve predictors versus the two or three for composition or abundance. Beyond day-of-year and plant composition, the simplified richness model included characteristics related to habitat diversity, weather, and location, suggesting a possible biogeographic effect on richness. As noted, canopy cover was a relatively more important predictor of butterfly community composition than of richness or abundance. For example, most observations of butterflies in the families Nymphalidae and Papilionidae were of individuals under canopy while Pieridae individuals were seen under canopy infrequently. A similar pattern of distribution was observed for Papilionidae and Pieridae in Sri Lanka [69].
Various other factors were important predictors of butterfly distribution. Temperature was a moderately important predictor of richness in regression tree models, suggesting a possible survey effect, such that more butterflies were observed between 22-32˚C than at higher or lower temperatures. Differences in the effective strip width of observations of different species also indicated how considering detectability affects descriptions of distribution. Specifically, we made up to nine-fold adjustments in calculations of species' density from raw counts due to differences among species in detectability. Comparing regression tree models, percentage of native plants present was a more important predictor of butterfly community composition than of abundance or richness. Local cover of tall vegetation (0.3-1 m tall) was most important as a determinant of butterfly abundance with abundance declining beyond low cover of tall vegetation. Finally, fire is a primary management tool for habitat restoration and maintenance in the mosaic landscape of habitats studied. Short-term (2-year) and long-term (15-year) fire frequencies were of moderate to low importance as predictors of butterfly community composition, richness, or abundance. In simplified models predicting composition, richness, or abundance, only long-term fire frequency was included as a predictor, and only for butterfly richness. However, the true effect of fire frequency on the butterfly community might be through the interaction of fire and plant community composition and canopy cover. For instance, in cork oak (Q. suber) stands in Portugal, removal of understory vegetation, which is often a goal of fire management, had a strong effect on butterfly richness and abundance that peaked within three years of mechanical removal of this vegetation [70] The butterfly community compositional differences among habitat types translated into significant differences in conservation value of the habitats for butterflies. Forests had the lowest conservation value for butterflies. Woodlands and scrublands had the highest conservation value, suggesting added conservation importance of landscapes with intermediate canopy cover and mid-successional stages-the conservation value index calculated as a measure of abundance of the most threatened butterflies in this study peaked at about 59% canopy cover. McCleery et al [18] examined how shifting woody vegetation density affected a trio of vertebrate groups (birds, bats, small mammals) in African savannas. These savannas were often undergoing changes that lead to canopy homogenization, either through loss of trees and shrubs or widespread increases in tree and shrub density. This replacement of the heterogeneous savanna canopy with more homogeneous low or high canopy conditions led to declines in diversity at the gradient extremes, illustrating a relationship between diversity of vegetation and structure and animal diversity [71,72]. However, in the Midwest USA, doubling of canopy cover over sixty years did not result in significant change in understory plant diversity but great species turnover [57]. Košulič et al. [21] examined the relationship between spiders and canopy cover in formerly coppiced woodlands, finding that a variety of spider community attributes-richness, functional diversity, activity density, and indicators of conservation concern-were related to canopy cover but not in a manner that desirable states of those attributes co-occurred at similar canopy covers. The greatest concentration of threatened and endangered spiders was found in range of canopies covers from 65-75%, indicative of woodlands. Studies of butterflies in a regenerating forest in Sri Lanka detailed that butterfly abundance and richness declined with increasing canopy cover with richness peaking near 20% cover [69]. Similarly, Ubach et al [73] documented replacement of grassland affiliated butterfly species with closed canopy species as forests expanded and replaced grasslands in northeast Spain with grazing and farming land abandonment. The native butterfly fauna was dominated by species associated with open habitats and the replacement of those open habitats with forests was associated with frequent local extinctions of butterfly species, most of which were of butterfly species that preferred open habitats. These are all examples of disparate animal taxa responding to heterogeneity in woody vegetation cover in ways that either abundance, richness, or conservation value were positively related to canopy heterogeneity at some spatial scale or that exhibited compositional change with changes in canopy cover but not changes in richness. In the Midwest USA butterfly community studied here, composition, richness, and abundance were all related to canopy cover but not all in the same way nor with the same strength of relationship. Thus, composition was more strongly predicted by canopy cover than was richness or abundance. Richness peaked at a lower canopy cover than did overall abundance. Conservation value peaked at intermediate canopy cover. Most of the butterfly species were not significantly associated with a single habitat type along the open to forest canopy gradient examined. This suggests that multiple habitat types might be used by the adults, perhaps for different resource needs, such as for oviposition, mating, and feeding [74]. Also, we examined habitat use by adults only and not by earlier life stages. While we documented a significant relationship between habitat diversity and species richness, conservation of habitat diversity might also be important for retention of individual species since different life stages might take place predominantly in different canopy defined habitats.
Compared to results from biodiversity studies of other taxa in our study area [12,13,28,75], the best predictors of different butterfly community attributes were not always the same as for those community attributes for other taxa but there were many similarities. For example, among the bird species that were significant habitat indicators most were associated either with the most open habitats, as we saw for butterflies, or the most forested habitats, which did not occur for butterflies. We found that bee abundance was negatively related to canopy cover, somewhat similar to butterflies, and positively to recent fire frequency, which was relatively an unimportant predictor of butterfly abundance [13]. Bee richness was positively related to plant richness while butterfly richness was related strongly to plant community composition but relatively weakly to plant community richness. Bee community composition was significantly related to plant richness and canopy cover, butterfly community composition was weakly predicted by plant richness and relatively strongly by canopy cover. Improving the status of savannas, and other habitats with intermediate canopy cover in the Midwest U.S., from their currently imperiled state [1,7,8], requires an understanding of how different taxa will respond to canopy manipulation. In studies of other taxa at the sites reported on in this study, habitats with intermediate canopy cover, such as savannas and woodlands had a high number of species relative to the canopy cover extremes, open and forest [11,13,28]. In this study, the species diversity patterns were somewhat different. Neither butterfly species richness nor species density peaked in savannas or woodlands. Species richness peaked in scrub habitats and species density in open habitats, suggesting that the number of butterfly species one would expect to find in a given plot on a given day was highest in open habitats but that the number of species present through time was greatest in scrub habitats. While savannas and woodlands in northwest Indiana might not be of the greatest value for maximizing the diversity of butterfly use, the value of woodlands might be found in higher use by butterfly species of conservation concern.
Conservation prioritization can be based on a variety of community attributes, such as diversity, abundance, or aiding threatened species [76,77]. For the Midwest USA butterfly community, these desirable conservation traits-maximizing diversity, increasing abundance, aiding the most threatened species through habitat characteristics-did not necessarily cooccur in one habitat type consistently. However, co-occurrence of higher diversity, overall abundance, and abundance of threatened species did increase as canopy cover increased to intermediate levels, ca. 67%, and then declined suggesting that habitats of intermediate canopy cover might be particularly effective for butterfly conservation in this region.