Effects of removing woody cover on long‐term population dynamics of a rare annual plant (Agalinis auriculata): A study comparing remnant prairie and oldfield habitats

Abstract Worldwide, grasslands are becoming shrublands/forests. In North America, eastern red cedar (Juniperus virginiana) often colonizes prairies. Habitat management can focus on woody removal, but we often lack long‐term data on whether removal leads to population recovery of herbaceous plants without seeding. We undertook a long‐term study (17 years) of numbers of the rare annual plant Agalinis auriculata in a gridwork of 100 m2 plots in adjacent prairie and oldfield sites in Kansas, USA. We collected data before and after removal of Juniperus virginiana at the prairie. Plant population sizes were highly variable at both sites and over time. High numbers of plants in a plot 1 year were often followed by low numbers the following year, suggesting negative density‐dependence. Plant numbers were lowest with extensive woody cover and with low precipitation. After woody plant removal, A. auriculata increased dramatically in abundance and occupancy in most years; increases were also seen at the oldfield, suggesting later survey years were overall more favorable. Synthesis and applications: Removal of woody plants led to increased numbers of a rare annual prairie plant, without seeding. Multiple years of data were essential for interpretation given extreme temporal variability in numbers. The largest prairie population was 7 years following tree removal, showing that positive effects of management can last this long. This species also fared well in oldfield habitat, suggesting restoration opportunities. Given that land managers are busy, time‐efficient field methods and data analysis approaches such as ours offer advantages. In addition to general linear models, we suggest Rank Occupancy‐Abundance Profiles (ROAPs), a simple‐to‐use data visualization and analysis method. Creation of ROAPs for sites before and after habitat management helps reveal the degree to which plant populations are responding to management with changes in local density, changes in occupancy, or both.


| INTRODUC TI ON
Temperate grasslands are conservation priorities; these once vast ecosystems have a long history of transformation to agriculture.
Grasslands around the world now face a new threat: conversion into shrublands or forests as a result of woody plant colonization.
Studies at the local spatial scale are particularly relevant to land managers who want to know how management decisions affect grassland populations and communities. For example, removal of eastern red cedar (Juniperus virginiana) can sometimes lead to recovery of herbaceous communities (Alford, Hellgren, Limb, & Engle, 2012;Limb, Engle, Alford, & Hellgren, 2014;Pierce & Reich, 2010).
These studies took a broad perspective, using many sites and taking data on percent cover or biomass of many species. Although invaluable from a community perspective, such data may lack adequate resolution for species-specific investigations. Yet ultimately the question of community transformation or recovery depends on the population ecology of individual species: that is, are increases or decreases of woody cover associated with predictable growth or decline of herbaceous plant populations? Dedicated population approaches are particularly needed for rare species that may not occur regularly in small sampling plots.
From a population ecology perspective, two complementary approaches are available: demographic studies and long-term surveys of plant numbers. The former consists of following plants to document survival and reproduction and developing demographic models (Caswell, 2001). For example, Andrieu et al. (2013) found that reproduction and asymptotic population growth rates of a rare perennial increased greatly after forest cutting and were comparable to open-habitat populations. Demographic models, however, often fail at forecasting dynamics of future years (Crone et al., 2013). This result is not surprising: environments change over time and these models are not designed to take into account all factors, including density-dependence.
In contrast to demographic models, long-term survey data document actual population dynamics-did a population decline in numbers with increased woody cover, and increase when woody plants were removed? Recording data in a continuous gridwork of plots allows analyses of both abundance and occupancy, and is important since plants in future years may occur in parts of the site where they are not currently found (Crawley, 1990). Survey data are most useful if the plant's life cycle is short relative to the data set length. It is thus not surprising that many studies of annual population dynamics have been done (Garcia de Leon, Freckleton, Lima, & Navarrete, 2014;Plaza, Navarrete, Lacasta, & Gonzalez-Andujar, 2012). Annual plants are challenging to study, however, because of yearly fluctuations in numbers, which can be due to processes originating within (endogenous) or outside (exogenous) of the population (Plaza et al., 2012). A common endogenous process is negative density-dependence (Garcia de Leon et al., 2014;Gonzalez-Andujar, Fernandez-Quintanilla, & Navarrete, 2006). Typical exogenous processes are temperature and precipitation (Garcia de Leon et al., 2014;Levine, Mceachern, & Cowan, 2011). Annual plants also often have seed banks, and germination of buried seed from past years may be important in population dynamics (Alexander, Pilson, Moody-Weis, & Slade, 2009;Salguero-Gomez, Siewert, Casper, & Tielborger, 2012).
In North American tallgrass prairies, the annual life form is uncommon: most plants are perennial. Our work examined the rare annual Agalinis auriculata in Kansas, USA in the context of woody colonization by eastern red cedar (Juniperus virginiana; Figure 1). This tree has increased across central North America (Meneguzzo & Liknes, 2015), primarily due to fire suppression (Ratajczak et al., 2014). As woodlands expand, grassland habitats decline: Briggs et al. (2002), for example, found fewer than five species/plot in highdensity cedar regions compared to 20-30 species in plots without cedar. Demographic effects are also apparent: Albrecht, Becknell, and Long (2016) noted that reproduction of the perennial Astragalus bibullatus declined in areas with cedar colonization.
Past work suggests that woody species have negative effects on A. auriculata. In Ohio, the largest populations were in areas with tree and shrub removal (Knoop, 1988). Based on a 4-year study in Illinois, Vitt, Havens, Kendall, and Knight (2009) concluded that populations should decline in the presence of woody brush (lambda = 0.81) and grow (lambda = 1.22) with brush removal. They cautioned that the F I G U R E 1 Flowering individual of Agalinis auriculata. Photograph by Craig Freeman positive effects of management on A. auriculata may be temporary, since competitors (often grasses) also increase after woody removal.
We began studying A. auriculata in 1996 at a remnant prairie, and expanded our work in 1997 to include an adjacent oldfield. In 2006, we removed J. virginiana and other woody species at the prairie with the goal of improving habitat for A. auriculata. Data collection continued through 2013, making it one of the longest data sets for terrestrial annuals outside of agriculture (Plaza et al., 2012). We first tested the hypothesis that woody plant presence alters abundance of A. auriculata. Second, we tested Vitt et al.'s (2009) hypothesis that positive effects of woody removal may be of short duration. Third, we addressed the degree to which A. auriculata is dependent on remnant prairie: were numbers comparable in the oldfield and prairie?
With both sites, we also explored how other factors were associated with plant numbers, ranging from precipitation to local variation in past plant numbers.
Our study was conducted at adjacent sites at the University In fall 2006 (after surveys, see below), all woody vegetation (except small shrubs) was cut in the prairie at ground level. Over 90% of the cut plants were J. virginiana. See Supporting information Appendix S1 for study site details.

| Plant surveys
In 1996, a 10 m by 10 m gridwork was overlaid on the prairie, producing

| Precipitation
Precipitation data (standard can measurements) were available as monthly summaries from a Field Station climate monitoring site 1 km southeast of the sites (https://biosurvey.ku.edu/field-station). We analyzed annual precipitation (cm) using a 12-month pe-

| Woody cover and eroded soil
We created a georeferenced shapefile of the grid using ArcMap10.3.1.
This grid was overlain on rectified aerial photographs to permit visual estimation of woody and eroded soil cover.

| Data sets
We created two types of data sets. First, to explore dynamics over For both sites, we excluded plots along the path or that had woody brush piles in some years. The average plant data set had 17 (prairie) or 15 (oldfield) data rows (1 row per year).
Second, we created a local data set to address plot-to-plot variation and to consider possible effects of neighboring plots. For the prairie, the local data set had 66 rows, corresponding to the 66 plots surveyed from 1996 to 2013 (excluding 2011). The oldfield local data set had 119 rows (119 plots surveyed from 2001 to 2010). These local data sets did not include plots along the path or those that had woody brush piles in some years. For each plot in both local data sets, we summed the number of plants in all 8 plots that surrounded it. Such "neighbor plants" could be a source of seed dispersal and thus contribute to numbers in a plot in the subsequent year. In counting neighbors, we included plots that were on the path or ones with woody brush piles since such plots had been surveyed for plants in all years except 2013. If a plot was on the gridwork edge, it had fewer than 8 neighbors; this was not problematic since the gridwork area covered the area where A. auriculata was abundant and we expected no (or only a few) plants outside the gridwork.

| Overview
With the exception of Rank Occupancy-Abundance Profile (ROAPS) analyses (see below), our goal was to explore whether variation in plant abundance within or across years was associated with exogenous factors (e.g., precipitation in a current [t] or previous year (t − 1); percent woody cover and percent eroded soil in the current year [t]) and endogenous factors (e.g., number of plants present in F I G U R E 2 Aerial photographs with grid (100 m 2 plots) before and after woody removal. The path (dashed line) separates the prairie (66 plots, north of path) and oldfield (119 plots, south of path). (a) 2002 color infrared image, before vegetation management; (b) 2006 color infrared image immediately after woody removal on the prairie; "X' marks locations where cut trees/ brush were piled; (c) 2015 color image after all management, including burning and mowing. Plots on the path or that had tree/ brush piles on them after woody removal were not used in most analyses; see Section 2.3. Light-colored patches of eroded soil are apparent in 2002 and 2015 prairie images.
the plot the previous year (t − 1) and number of neighbor plants in both the current year t and previous year t − 1). Details are described below; in all cases, we first fit a global model and compared the fit of models with fewer parameters to the global model. We used Akaike's information criterion (AIC) values to select the most parsimonious models from the set of explored models. If a variable was kept in the model, the most important feature is the sign of the coefficient (i.e., if positive, that variable was associated with increased plant numbers). Following Burnham and Anderson (1998), the best model had the lowest AIC score and we considered models with a ΔAIC > 2 to have more support than other competing models and used parameter estimates derived from that model. We used SAS (version 9.4).

| Average number of plants/plot across years
We used multiple linear regression to determine what variables best predicted the average number of plants/plot. We used the average plant data sets in separate analyses for the prairie (1996-2013, excluding 2011) and oldfield (1997-2013, excluding 2011, 2012). Variables in the initial model for each site were precipitation t , precipitation t-1 , number of plants t-1 , number of neighbors t , and number of neighbors t−1 . Average number of plants per plot was log-transformed to satisfy the assumption of normality in linear regression.

| Number of plants/plot: within years (prairie)
We also determined which variables best predicted plant numbers within a year. We used log-linear models in single year analyses with only the local prairie data set. We analyzed years before (1997,1999,2001) and after (2008, 2009, 2010, 2012, and 2013) woody removal (only years where at least one plot had>5 plants).
Variables included in initial models were percent woody cover t , percent eroded soil t , number of plants t-1 , number of neighbors t , and number of neighbors t−1 . We compared fit of zero-inflated models with Poisson models (not zero-inflated). Models for 2009 and later were best described by a log-linear model with a Poisson distribution. Models before 2009 were best fit with a zero-inflated Poisson.

| Comparison: 2001-2006 versus 2007-2010 (plants/plot)
Using the local data sets, we ran four models: both the prairie and oldfield sites for 2001-2006 (before woody removal in the prairie) and for 2007-2010 (after woody removal in the prairie).
Woody cover did not change over time at the oldfield, but we compared both sites for the same periods to evaluate if factors other than woody removal were changing plant populations.
Variables included in the initial model were precipitation t , precipitation t-1 , number of plants t-1 , number of neighbors t , and number of neighbors t-1 . We used a zero-inflated model with mixed effects to fit regressions (plot was a fixed effect; year was a random effect). We compared fit among possible distributions (negative binomial or Poisson) by comparing the predicted number of zeros with the observed. All models were best fit by a zero-inflated Poisson model.

| Comparision: 2001-2006 versus 2007-2010 (abundance and occupancy)
We used the local data sets (prairie and oldfield) and Rank Occupancy-Abundance Profile (ROAP) analyses (Collins, Holt, & Foster, 2009) to explore temporal changes in local abundance (number of plants/plot) and occupancy (proportion of plots occupied). In a ROAP plot, the y-axis displays local abundance in each plot, and the maximum value represents the plot with the highest abundance. Plots are ranked in order of their abundance along the x-axis, with the highest abundance plot being ranked first ("1").
When the rank is divided by the number of plots surveyed, the maximum x-value displays the proportion of plots occupied in the landscape (i.e., occupancy). A ROAP plot for a year allows one to quickly visualize the maximum abundance (y-intercept), the range of abundances across the plots (other y values), and the occupancy (x-intercept). Total abundance is the area under the ROAP; the area between two ROAPs displayed on the same graph reflects magnitude of change in overall abundance, accounting for both spatial expansion or retraction, as well as local plot densities. See Supporting information Appendix S2 for example data sets and ROAPs.
We constructed ROAPs for both prairie and oldfield for each year. To reduce the number of statistical comparisons and to focus on our research questions, we also constructed average ROAPs for two times: before woody removal (T 1 ), and after woody removal (T 2 ). We made two comparisons: (a) using the same years where data were available in both prairie and oldfield (before, 2001-2006 vs. after, 2007-2010) and (b) using all data (prairie only; before, 1996-2006 vs. after 2007-2013, excluding 2011).
The former allowed us to explore if factors other than woody plants were altering population dynamics. Averages were calculated for each rank over the appropriate time periods on data sets where each year had been sorted by plot abundance. In plots of average ROAPs, the y-axis reflects the average local density in occupied plots. Because density was averaged across multiple years, the maximum x-value is the maximum occupancy achieved over the years under investigation (Supporting information Appendix S2).
We used randomization tests to determine whether total abundance shifted over time at each site. Specifically, we randomly assigned local densities to a time period, calculated the area between the two ROAPs (D*; Collins et al., 2009) 999 times, then compared the empirical D* to this distribution. We considered the change in abundance statistically significant if the empirical value for D* was among the 50 most extreme values.

| Overview
Numbers of A. auriculata plants varied widely at both sites, both temporally ( Figure 3) and spatially (Supporting information Figure   S1). Plots that had most plants in 1 year were not always the plots that had most plants in other years (Supporting information Figure   S1). Some years had relatively few plants at both sites (2010) Figure   S1). Except for 2007, these years also had high woody cover (Figures 2a, 3). The average amount of woody cover at the prairie

| Average number of plants/plot: across years
Average number of plants/plot was marginally higher in the oldfield than the prairie for all years where data were available (paired t test, t 14 = 1.94, p = 0.072); there was no difference after woody cover was removed (although sample size was low; paired t test, t 4 = 1.13, p = 0.320). For the prairie, the average number of plants/plot was best predicted by average percent woody cover. The model with the lowest AIC only included woody cover, with higher plant number associated with lower woody cover: (log N t = 0.90-0.03 average percent woody cover; Figure 3). In

| Number of plants/plot: within years (prairie)
Prior to woody removal, greater woody cover was associated with fewer plants (2 of 3 years) and a higher percentage of eroded soil was associated with fewer plants (1 of 3 years; Table 1A). The models consistently showed a positive association between numbers in the previous year and numbers in the current year. Usually, the number of neighbors (previous or current year) were negatively associated with plant numbers.
After woody removal, increased woody cover was always negatively associated with plant numbers, though not included in all models ( were positively or negatively associated with number of plants. Table S1 for AIC values for the best and second best models for models in Table 1A,B, as well as ΔAIC values.  Table S2 for AIC values for the best and second best models for models in Table 2, as well as ΔAIC values.

| Abundance and occupancy
There was considerable year-to-year variability in abundance and occupancy at the prairie both before ( Analyses performed with the prairie data for all years (before: 1996-2006, after: 2007-2013, excluding 2011) showed similar increases in abundance and occupancy (Supporting information Figure   S2, D* = 889.53, p < 0.001).

| Factors affect plant population numbers
Our long-term study of the rare annual Agalinis auriculata revealed that woody colonization and removal were key factors affecting numbers. Average woody cover was the only predictor of average TA B L E 1 Parameter estimates (and SE) from best-fit models describing how numbers of Agalinis auriculata in 66 prairie plots (each 100 m 2 ) depend on percent woody cover, percent eroded soil, numbers of plants in the previous year, neighbors ( number of plants per prairie plot over a 17-year period, with plant numbers declining as woody cover increased. Prairie plots with higher woody cover had fewer A. auriculata plants in most years. We also note that both prairie and oldfield had more similar population models after woody removal than before, suggesting that dynamics at the two sites converged after woody removal. Finally, although plot occupancy and plant abundance were higher after than before woody removal at both sites, the magnitude of increase over time was much greater for the prairie (in part because it was so rare prior to tree removal). Overall, these results suggest that the A. auriculata We lack data on what factors affected woody colonization at the site, as well as the mechanism for decline in A. auriculata when J. virginiana became common. Agalinis auriculata was likely shaded out by trees, but water availability or soil properties could also have been altered. Given that A. auriculata is a hemiparasite, another possibility is that the high tree cover reduced host numbers, although probable hosts (Asteraceae; Molano-Flores et al., 2003) were present at both sites and over time.
Other exogenous factors also were associated with plant numbers, especially precipitation. Not surprisingly, high precipitation years had high plant numbers. We do not know in detail how precipitation affected plant survival or reproduction; in some annuals, for example, the timing of rain was more important than total annual rainfall for population dynamics (Levine et al., 2011). Different factors may also interact: Figure 3 shows that the years immediately following woody removal also had high rainfall, likely facilitating F I G U R E 5 Rank occupancyabundance profiles (ROAPs) for Agalinis auriculata. ROAPs were constructed for (a) prairie (individual years), (b) prairie (averaged across years for two time periods: before woody removal in the prairie (red, 2001-2006) and after woody removal in the prairie (blue, 2007-2010), (c) oldfield (individual years), and (d) oldfield (averaged across years for two time periods: before woody removal in the prairie (red, 2001-2006) and after woody removal in the prairie (blue, 2007-2010). For plots of individual years, local abundance was measured as numbers of plants in a 100 m 2 plot (y-axis) and the x-axis refers to the relative rank (i.e., a plot with the highest abundance has the lowest relative rank Endogenous processes, including negative density-dependence, are often important in plant population dynamics (Crawley, 1990;Garcia de Leon et al., 2014;Gonzalez-Andujar et al., 2006;Plaza et al., 2012). Consistent with this concept, plots with large numbers of plants in 1 year had smaller numbers of plants in the next year and vice versa. Using the average plant data set, the oldfield model has a positive term for precipitation and a negative coefficient for previous year numbers. Perhaps high water availability contributed to an increased population size, but plants at high density had reduced per-capita seed production. Multiple mechanisms could be at work: intraspecific competition, disease spread at high densities, or concentrated herbivore feeding in high-density plots are possibilities.
We did not directly test for density-dependence because our local data set included plots without plants (zero density; no possibility for density-dependent processes). Interestingly, the coefficient for previous year numbers was often positive in these models, perhaps suggesting that the presence of plants in past years was indicative of conducive locations for plant growth in future years.
These positive coefficients may also be a byproduct of using zero-inflated binomial distributions, which in effect may reduce the role of the zeroes in the analyses.
Plots with zero plants in 1 year sometimes had many plants present in subsequent years. These recruits could be from seed dispersal, from germination from dormant seeds (Baskin et al., 1991)  year before occurred in studies of annual sunflower, another seed bank species (Alexander et al., 2009).
Although we expect that seed banks, as opposed to seed dispersal, are the most likely explanation for transitions from zero plants/ plot to 100's of plants/plot, as noted above, seed dispersal is also likely ecologically important. We were surprised to find A. auriculata at the oldfield site in 1997 given that the area had been tilled in 1995 (Supporting information Appendix S2); presumably these plants are the result of dispersal from the nearby prairie. We do not know the degree to which the numbers of plants at the oldfield in later years were derived from dispersal from the prairie or a result of seed production at the oldfield site.
On a small spatial scale, we also explored the role of seed dispersal when we added the number of neighbors to models to include some spatial structure because we expected that neighbors in a previous year could disperse seed into a plot and thus be positively

| Rarity, restoration, and roaps
In addition to studying A. auriculata in a remnant prairie, our work in a post-agricultural field provides data relevant to restoration.
Specifically, although known as a rare plant, A. auriculata colonized an oldfield and persisted with often large population sizes. This species thus has a tolerance of disturbed sites, and potentially could be introduced into prairie restoration plantings, especially since Asteraceae (hosts for this hemiparasite) are common in prairie seed mixes. As noted previously, we found that seeds can at least sometimes disperse large distances (i.e., meters as opposed to centimeters), suggesting a "spillover" effect where a desirable species disperses from a remnant site (prairie) into non-target habitat (oldfield; Brudvig, Damschen, Tewksbury, & Haddad, 2009).
Like most ecological studies, our data are site and speciesspecific. However, our work suggests three general lessons for management of herbaceous grassland species where woody colonization is occurring. First, we show that improvement in population trends for a rare species is possible with tree cutting and without reseeding: too often managers may not take any action because they feel any efforts are hopeless. Second, as perhaps obvious, long-term records are essential, given that current threats (woody colonization) and management responses (woody removal) occur over many years. Our results (Figure 3) reveal the tremendous fluctuation in numbers that can occur across years for annual species, illustrating that results from any 2-3 year period could lead to misinterpretations. Third, we note that by taking data with a relatively simple method (number of plants in large plots), we could extensively sample two adjacent habitats (prairie and oldfield). If we had only taken data in a subset of the area using small plots, we could easily have missed major trends, especially because plants "moved around" both habitats over the years (Supporting information Figure S1). Haddad et al. (2008) suggest that optimal monitoring schemes for rare butterflies require consideration of the need for information, statistical performance, and the costs of different approaches. These same factors must be considered with plants. With annuals, count data, such as presented here, have logistical advantages. Annual surveys required only 2-3 days with an established gridwork. Of course, combining demographic and survey approaches would be ideal and allow greater understanding of the mechanisms leading to changes in numbers of plants. For example, plant size and reproduction data would have allowed us to compare our work to Vitt et al. (2009) who noted that the proportion of large individuals (who produce the most seed) increased with woody removal. Regardless of monitoring approach, it is obvious that one needs consistency in data collection methods across years and careful record-keeping. In our case, work had been largely done by the same people for the 17 year period, but that is rarely the case and many long-term studies are doomed by poor data quality.
ROAPS are another tool to consider in optimal monitoring schemes. Although not spatially explicit, ROAPS provides insights into spatio-temporal shifts in abundance that may not otherwise be apparent. It is also important to emphasize that many land managers have limited time, and monitoring schemes does not always require the most advanced tools. Simple monitoring of local plant densities and spatial spread, as done in this study, can be used to plot ROAPs (single year or averaged over years, Supporting information Appendix S2) to compare abundances over time, between landscape types, or before/after management actions. Such visual tools allow land managers to discern whether overall abundance is shifting due to changes in local density, spatial extent, or both, and thus more readily use past data to guide future management decisions. For instance, ROAPs could be useful for invasive plant research: management responses would likely differ if occupancy increased more rapidly over time than did local abundance, as opposed to situations where the plant was highly abundant locally but not spreading.

ACK N OWLED G M ENTS
We particularly thank Galen Pittman for his careful and conscientious data collection and work on habitat restoration; many students also assisted. Jorgina Ross helped interpret photographs. We appreciate discussions with Jim Bever, Craig Freeman, Dan Hirmas, Kelly Kindscher, and Ruth Shaw.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
WDK, HMA, and CDC conceived ideas and designed methodology; VBS, WDK, HMA, and DAC collected data; AWR, CDC, HMA, WDK, DAC, and LDC analyzed data with AWR leading on general linear model analyses and CDC leading on ROAPs analyses; HMA and CDC led writing of the manuscript. All authors contributed critically to drafts and gave final approval for publication.

DATA ACCE SS I B I LIT Y
Data available from the Dryad Digital Repository: https://doi. org/10.5061/dryad.801r9q8.