Study Area

A P. praeclara population exists at Pipestone National Monument in southwestern Minnesota (44° 0.799 N, 96° 19.518 W). The population occurs in mesic prairie habitat that was part of the original 46 hectares designated as a national monument in 1937, which had escaped cultivation. The soil type consists of Ihlen and Ihlen-Rock outcrop complex (Becker et al. 1986). Ihlen soils are well drained and consist of 23–40 inches of loess over Sioux quartzite on relatively level surfaces. The Ihlen-Rock complex includes rock outcrops within thinner Ihlen soils. Shallow drainages are found in both soil types. Wildfires occurred until the early 1940s, when relatively effective fire suppression began. Resource managers initiated prescribed burning in the orchid habitat in 1973.

Sampling Methodology

Surveys were conducted for P. praeclara annually from 1993 to 2012. Monitoring took place during peak flowering, which always occurred in July, with the exception of 2012, in which peak flowering occurred during the last week of June. In 1993 and 1994, a haphazard search of the area was conducted by two observers. Beginning in 1995, a systematic search of the area was conducted by a team of four to eight observers walking in tandem 2 m apart (Willson 2000; Young et al. 2007). Locations of flowering orchids were recorded with a GPS unit. Accuracy of GPS positions was within 5 m through 2006, and within 2 m for 2007 and later years. Plant height and the number of flowers and buds were also recorded. Given the more thorough nature of the surveys beginning in 1995, the number of flowering orchids for 1993 and 1994 could be underestimates, and plant height and number of flowers/buds could also be biased. Thus, we include these data for descriptive purposes, but for modeling effects of precipitation we include only years in which the same monitoring methodology was used (as did Willson et al. 2006). Prescribed bums were conducted in four years of this study: 1994, 1997, 2002, and 2009, always in May. The entire orchid habitat was burned as uniformly as possible.

Statistical Analyses

We evaluated the effect of precipitation on number of flowering plants during six phenological stages as identified by Wolken (1995) (Table 1). We first tested the model obtained by Willson et al. (2006) (referred to as the Willson et al. model in the results): (Number of flowering plants = 196.73 + 7.28 precipitation during senescence/bud development - 9.30 precipitation during dormancy), which explained 77% of the variability in flowering plant number for the years used to construct the model. We compared the observed number of flowering plants to the number predicted by the model in each of the eight years following publication of the model, and calculated the average absolute difference between observed numbers and the model predictions. We also determined for each year, as a null model of no precipitation effect, the difference between the model prediction and the average number of flowering orchids over the same eight years on which the model was based. We calculated the average absolute difference between observed numbers and this longterm average.

Next, we re-evaluated the relationship between number of flowering plants and precipitation in the specified phenological stages by employing an information-theoretic approach. In this method, inference is based on Kullback-Liebler (K-L) information (Kullback and Liebler 1951). Akaike's Information Criterion (AIC) (Akaike 1973) was used as an estimator of relative K-L information (Burnham and Anderson 2002). We conducted all subsets regression, and used AIC to evaluate the regression models. Because the sample size was small, we calculated AIC adjusted for small sample size (AIC_c) as:

Table 1.

Six phenological stages of P. praeclara as identified by Wolken (1995). Time period refers to the dates over which precipitation was accumulated for each stage.

where K is the total number of parameters in the model (including the intercept and residual variance) and n is the sample size. Because AIC_c is on a relative scale, differences in AIC_c values were calculated as:

for all candidate models, where i indicates the i^th model. The best approximating model was the one with the lowest AIC_c score. Models with ΔAIC_c ≤ 2 have “substantial” support, however, and would serve nearly equally well in approximating the information in the data set (Burnham and Anderson 2002). We report all models with ΔAIC_c ≤ 2.

A total of 18 years of data obtained with a consistent methodology was available (1995–2012), although prescribed fires occurred in three years (1997, 2002, and 2009) (Appendix). Because fires may influence flowering of the orchid, an analysis was conducted excluding the burn years (n = 15 years), and a separate analysis was conducted including all years of data (n = 18).

We conducted parallel analyses using plant height and number of flowers/buds as response variables in place of number of flowering plants, to evaluate whether there was any evidence that precipitation in the six phenological stages affected these variables. As in the analyses focusing on number of flowering orchids, analyses for plant height and number of flowers/buds were conducted for all years of data (1995 to 2012) and excluding bum years. Because no orchids were present in 1998, this year was excluded from all analyses of plant height and number of flowers/buds. Data from 1993 and 1994 were not included; because of the haphazard nature of the sampling in those years, shorter orchids, or those with fewer flowers, may have been missed.

Durbin-Watson tests were used to evaluate the linear models with substantial support, to determine whether the error terms were serially correlated. Positively autocorrelated error terms often result from a model that is missing key variables (Kutner et al. 2005).

We evaluated whether the prescribed bums affected number of flowering orchids, plant height, or number of flowers/buds with t-tests and one-way analysis of covariance (ANCOVA). We evaluated potential effects of bums both in the year of the bum and the following year, in separate analyses, since flower bud initiation occurs in the year before flowering (Bowles 1983), and removal of litter from fire may affect orchids in the year after the bum (Pleasants 1994).

We obtained daily precipitation records for the Pipestone weather station (216565) from the National Oceanic and Atmospheric Administration (NOAA) Global Historical Climatology Network (GHCH) ( http://gis.ncdc.noaa.gov/map/viewer/). This weather station was located ∼600 m from the orchid population. A small fraction of the precipitation data was missing—0.6% of the total records—and in these cases, the average of three surrounding weather stations was substituted. The three weather stations used were Tyler (218429; ∼20 miles to the northeast), Flandreau (392984; ∼13 miles to the west), and Lake Wilson (214534; ∼18 miles to the east). Fifty-five percent of the missing values were zeros.

All subsets regression and AIC evaluation were conducted using SAS ® software (SAS Institute Inc. 2008). Inspection of diagnostic plots revealed no major departures from regression assumptions such as normality or constancy of error variance. Durbin-Watson tests, t-tests, and AVOVAs were done with IBM SPSS Statistics software (IBM Corp. 2011). The sequential

Bonferroni method (Rice 1989) was used to correct for multiple comparisons when necessary.

Demographics

The spatial distribution of the orchid was generally clustered around two low, wet areas (Figure 1). The number of flowering orchids was variable, but generally <150 were observed in the first 18 years of the study, with >200 individuals found in only one of those years (2003) (Figure 2A). Numbers were often quite low, with no orchids being found in 1998, and ≤3 in 1997 and 2007. In contrast, the number of flowering orchids in both of the final two years of the study was >600. The mean number of flowers/buds varied between 5 and 10 across the entire data record, but revealed a relatively steady decline over the first 17 years of the study, before rebounding in the last three years (Figure 2B). Mean plant height was greatest in the first four years of the study, and highly variable (albeit lower) afterwards (Figure 2C).

The number of flowering orchids was not significantly correlated with either mean plant height or mean number of flowers/ buds. Mean plant height and mean number of flowers/buds, however, were positively correlated (r = 0.49; P = 0.032), although this was not significant if corrected for multiple comparisons by the sequential Bonferroni method (n = 3 comparisons).

Effects of Precipitation

The Willson et al. model yielded predictions that were usually not close to the observed numbers of flowering plants (Table 2). Interestingly, burn years were excluded from the data used to construct the Willson et al. model, yet the most accurate prediction occurred in 2009, the only burn year among these comparisons. The average number of flowering plants in the data set used in the Willson et al. model was 87. If one were to use this simple average as a predictor for the number of flowering plants in years succeeding the model (2005–2012), the average absolute difference between the observed number and the average would be 187. The average absolute difference between the observed number and the number predicted by the Willson et al. model for this time period is 194. Thus, over this time period, simply using the historical average would be a slightly better predictor of orchid numbers than use of the model.

Figure 1.

Map panels showing distribution of flowering orchids in each of 20 years (1993–2012) at Pipestone National Monument, Minnesota. Accuracy of GPS positions was within 5 m through 2006, and within 2 m for 2007 and later years.

When the larger data set was analyzed by an information theoretic approach, and burn years excluded (n = 15 years), the best approximating model based on AIC_c scores included two categories: (1) mature growth in the previous year (the year before flowering); and (2) emergence (Table 3). Two other models, however, were essentially as good. These three models always contained some subset of three categories: (1) mature growth in the previous year; (2) postsenescence; and (3) emergence (or all three). The model containing the predictor variables senescence/bud development and dormancy (the two variables in the Willson et al. (2006) model) had ΔAIC_c = 16.42. A model with ΔAIC_c > 10 has “essentially no empirical support” (Burnham and Anderson 2002).

If burn years are included (n = 18 years), the same best approximating model was obtained (with somewhat different coefficients) (Table 3). Only one other model (with all three of the categories) was essentially as good, however. The relatively small differences in the coefficients are likely to result more from model selection uncertainty due to small sample sizes rather than an effect of fire. Again, the model containing the predictor variables senescence/bud development and dormancy performed poorly, with ΔAIC_c = 19.26.

Figure 2.

(A) Number of flowering plants from 1993 to 2012 at Pipestone National Monument, Minnesota. (B) Mean number of flowers and buds (Error bars = 1 SE). (C) Mean plant height in cm (Error bars = 1 SE). In all panels, open circles indicate years in which a prescribed burn occurred. Closed circles indicate years without fire.

All Durbin-Watson statistics were in the indeterminate range, yielding inconclusive test results (D ranged from 0.95 to 1.30). In principle, this is due to lack of enough cases to obtain a definitive result.

Analyses of models with plant height and number of flowers/buds as response variables revealed that relatively little of the observed variability in either response variable was explained by variability in predictor variables in the most parsimonious models. Models with ΔAIC_c ≤ 2 yielded coefficients of determination of ≤0.38 for number of flowers/buds and ≤0.24 for height in years with no burns (n = 14); maximum values of R² for all years (n = 17) were ≤0.14 for number of flowers/buds and ≤0.42 for height. We do not present the models here since all were relatively poor.

Effects of Fire

There was a marginally significant (P = 0.088) effect of fire on the number of flowering orchids, which was lower in burn years than nonburn years (43.3 ± 54.9 vs. 152.9 ± 216.7, mean ± SD). Because precipitation was observed to have important effects on the number of flowering orchids in the regression models, precipitation values during the three phenological stages found to be important (mature growth in the previous year, emergence, and post-senescence) were included as covariates in an ANCOVA. When these variables were taken into account, there was no effect of fire (F = 0.09, P = 0.77, df = 1).

The number of flowers/buds was not different in burn years vs. nonburn years (P » 0.05). Plant height was shorter in years with a prescribed burn (42.9 ± 9.6 vs. 52.7 ± 5.8 (mean ± SD); t = -2.63, df = 17, P = 0.018). If only years with a consistent, systematic monitoring are analyzed (all years except 1993 and 1994), this difference in plant height is very robust (38.1 ± 0.6 vs. 52.2 ± 5.7 (mean ± SD); t = -4.20, df = 15, P = 0.001).

Table 2.

Observed number of flowering orchids, number of flowering orchids predicted by the Willson et al. (2006) model, the absolute difference between observed and predicted numbers, and the absolute difference between observed numbers and the average number of flowering orchids (87) over the period covered by the Willson et al. (2006) model, at Pipestone National Monument, Minnesota. A burn occurred in 2009.

The three years with the shortest plants (mean height <40 cm) were all bum years (1997, 2002, and 2009); however, they also happened to be the three years with the least spring precipitation (Figure 3). In fact, there was a significant positive association between plant height and precipitation during the emergence phenological stage (r = 0.47, P = 0.043). If plant height as a function of prescribed bums is compared using an ANCOVA with precipitation during the emergence stage as a covariate, the effect of fire is no longer significant (F = 2.97, P = 0.104, df = 1).

Table 3.

The best approximating model (lowest AIC_c score) and competing models having “substantial” support (ΔAIC_c ≤ 2) for the P. praeclara population at Pipestone National Monument, Minnesota. MLL is the maximized log-likelihood. K is the total number of parameters in the model (including the intercept and residual variance). The response variable for all models is number of flowering orchids.

When similar comparisons were made in the year after fire (i.e., comparing bum years + 1 to all other years), there were no significant differences in either number of flowering orchids, number of flowers/buds, or plant height (all P > 0.05).

Effects of Precipitation

Based on an examination of P. praeclara root systems, Wolken (1995) proposed a phenological system of above- and belowground plant development. Our model results indicate the relative importance of precipitation in the various life stages of the previous and current growing season. The importance of precipitation during the previous growing season (June–July) suggests the existence of a lag effect. Some threshold of precipitation, or more likely soil moisture, may be required to promote robust growth of the perennating bud and supporting root system that leads to flowering during the subsequent season. As a more immediate effect, there may exist a precipitation threshold during emergence (April–May) that is required to further support plant development leading to flowering. The periods of mature growth and emergence correspond with periods of greatest photosynthetic activity in P. praeclara. This putative mechanism may also explain variation in P. praeclara population size (Pleasants 1995) and demographic transitions between flowering, nonflowering, and dormant states (Sieg and King 1995). While our modeling cannot ascertain the relative strengths or interactions of precipitation during these periods, we have detected a signal indicating the relative importance of precipitation that presumably corresponds with distinct phenological stages in P. praeclara.

Figure 3.

Mean P. praeclara height as a function of precipitation during the emergence stage from 1993 to 2012 at Pipestone National Monument, Minnesota. (Error bars = 1 SE).

Why did the models obtained using the information theoretic approach differ so greatly from the model obtained by Willson et al. using stepwise multiple linear regression? Numerous issues exist with stepwise model selection procedures, including the problem of multiple hypothesis testing, bias in parameter estimation, inconsistencies among model selection algorithms, and an inappropriate focus on a single best model (Anderson and Burnham 2002; Austin and Tu 2004; Whittingham et al. 2006; Wiegand 2010; Hegyi and Garamszegi 2011).

The most important factor, however, probably has to do with the smaller sample size (n = 8) evaluated by Willson et al. Different models may be determined to be “best” for different data sets (from the same population), and this source of model selection uncertainty becomes greater as sample size becomes smaller. The Willson et al. model had essentially no empirical support compared to the best models determined with larger sample sizes (n = 15 and 18). While our sample sizes were larger than those available to Willson et al., they are still relatively small. Instead of constructing a model based on a subset of the data and then testing it on another subset, we used all available data in model selection, in an attempt to obtain the most robust inference possible. Thus, although the models generated may fit these data fairly well, there is still the chance that they are not representative of the overall population.

The available data did not support selection of a single best model. Other methods, such as stepwise regression, are guaranteed to produce a single best model, even though there may exist substantial model selection uncertainty. The inability to identify a single best model is not uncommon when using an information theoretic approach, and reflects the natural complexity of the system, rather than a shortcoming of the method. Robust inferences can still be made from a set of competing models with substantial support, that is, multimodel inference (Burnham and Anderson 2002). More formal procedures include model averaging, although in this case model averaging would probably not be very useful. The inability to rule out positive autocorrelation in the error terms suggests that remedial measures for autocorrelation (see Kutner et al. 2005) should be explored and transformed models evaluated before attempting to make any inferences from parameter coefficients. Given the caveats above, identifying the important phenological stages, rather than specifying exact parameter coefficients, represents the most important inference from the analyses presented here.

Because of the relatively small sample size, and the fact that we evaluated all possible models, the results of this study should not be considered as confirmatory. Other environmental variables are likely to be important and, like most studies, we consider only simple linear models, although the true relationship may be curvilinear or contain higher-order interactions. The inconclusive results of the Durbin-Watson tests do not preclude the existence of other key variables. Thus, the results of the study should be taken more as an early step toward evaluating the relationship of flowering orchids with environmental variables such as rainfall during different phenological stages, and demonstrating that an earlier published model is not useful for prediction.

A limitation inherent in the information theoretic approach is that this method will only select the best model(s) from a larger set of candidate models. If the candidate models are all poor, then the selected model(s) will also be poor. In the analyses of number of flowering orchids, coefficients of determination (R²) were always ≥0.7, indicating a relatively large amount of the variability in the number of flowering orchids could be explained by variability in a few precipitation categories. In contrast, in the analyses of plant height and number of flowers/buds, coefficients of determination were always ≤0.42, indicating that the majority of the variability in these two variables was due to other factors outside the scope of these analyses. Finally, our analyses considered only one P. praeclara population. Thus, these results may not be directly generalizable to other populations in other areas.

Effects of Fire

Anecdotal observations in the literature suggest that fire may increase flowering in P. praeclara (Sheviak 1974; Bowles 1983; Currier 1984). The effect of fire was difficult to evaluate, since only four prescribed bums occurred over the course of the data record. Moreover, fire in this case is a complex variable, as both the timing and intensity of the bums varied. No consistent increases or decreases were observed in orchid numbers following burns. Because the best models changed little when fire years were included, it can be inferred from the regression analyses that fire, as applied during the month of May, did not greatly affect how the number of flowering orchids responded to precipitation in the life stages evaluated. The ANCOVA analyses provided complementary results, finding no significant effect of fire on the number of flowering orchids after taking into account the variability in precipitation.

Shorter plants were observed in years in which prescribed fire occurred, although most of the burns occurred in years with low spring precipitation. Pleasants (1994) conducted a controlled burning experiment and hypothesized that growth in P. praeclara may be affected by an interaction of fire and soil moisture, with fire in dry years suppressing growth rates, and fire in wet years increasing growth rates. His study was limited, however, in that it only included bums in two years. While our results are supportive of this hypothesis, a more parsimonious interpretation of our data is that precipitation was important in determining plant height, but fire was not, as the effect of fire disappeared when variability in precipitation was taken into account. We also found no evidence for a one year lag in the effect of fire.

Thus, overall, our data support important ecological effects of precipitation but no effects of fire applied in May on orchid growth or reproduction. An important caveat is that if fire intensity, fire timing, and precipitation do interact (and these factors may exhibit considerable variability), then a much larger sample size would be required to evaluate such interactions.

Land managers face a dilemma when using prescribed fire to maintain habitat occupied by P. praeclara, while preventing direct injury to orchids. The prescribed fires during May evaluated in this study were timed to reduce the abundance of the exotic cool season grass Bromus inermis Leyss (smooth brome; Willson and Stubbendieck 1997). Late spring fires may damage P. praeclara individuals likely to flower, however, and fall bums provide an alternative strategy for reducing smooth brome (Biederman et al. 2014). While no effect of fire was observed on the orchid populations of our study, relatively few bums occurrred over the course of the data record, and no control plots existed. The timing of emergence of P. praeclara is variable, and spring rains often make burning difficult. Thus, bums in October and November, when the orchid is dormant, may provide better opportunities to introduce fire without impact to P. praeclara, although this may increase the probabililty of smooth brome expansion.