Assessing the hierarchy of long-term environmental controls on diatom communities of Yellowstone National Park using lacustrine sediment records

Abstract

Chraïbi VLS, Fritz SC. 2020. Assessing the hierarchy of long-term environmental controls on diatom communities of Yellowstone National Park using lacustrine sediment records. Lake Reserv Manage. 36:278–296.

An ecosystem’s ability to maintain structure and function following disturbance, defined as resilience, is influenced by a hierarchy of environmental controls, including climate, surface cover, and ecological relationships that shape biological community composition and productivity. This study examined lacustrine sediment records of naturally fishless lakes in Yellowstone National Park to reconstruct the response of aquatic communities to climate and trophic cascades from fish stocking. Sediment records of diatom algae did not exhibit a distinct response to fish stocking in terms of assemblage or algal productivity. Instead, 3 of 4 lakes underwent a shift to dominance by benthic diatom species from 1935 to 1950, which suggests lower lake levels resulting from warmer, drier climatic conditions. The lake that did not undergo such a shift is fed by groundwater rather than snowmelt, suggesting a buffering effect by water source. Dissimilarity analysis shows that the diatom assemblages in all 4 lakes have become increasingly dissimilar since circa 1955, suggesting that communities have not yet stabilized from the first-order influence of climate. Thus, climate likely had a more prominent influence on diatom community structure than did manipulation of the fish community. Understanding the relative importance and interplay among abiotic and biotic stressors and the resultant resilience of an ecosystem provides implications for the adaptive management of lakes.

Keywords:

• Climate change
• diatoms
• fish stocking
• paleolimnology
• resilience
• trophic cascade
• Yellowstone National Park

Holling (1973) first proposed resilience theory as a framework for understanding how ecosystems respond to environmental stress. Although the definition of resilience varies, here we define ecological resilience as the ability of an ecosystem to withstand disturbance and to resist transitioning into an alternate stable state with different ecological structure and function (Gunderson 2000). Ecological stress affects a system over multiple temporal and spatial scales; hence, a major regime shift can result either from a single large disturbance or from the cumulative effects of several smaller scale changes (Gunderson 2000). An addition to resilience theory, panarchy, characterizes the interconnected, multiple-scale processes that determine the fundamental resilience of a system (Walker et al. 2004).

An obstacle to the study of resilience is the difficulty of conducting research over the long temporal scales needed to capture empirical evidence of the hierarchy of environmental influences that characterize panarchy, which is necessary to apply resilience theory to adaptive management of natural systems (Gunderson 2000, Walker et al. 2004, Fritz and Anderson 2013). Paleolimnological research has the potential to provide empirical evidence of resilience (e.g., Spanbauer et al. 2014) by extending the temporal resolution of ecosystem dynamics to generate a holistic picture of multiple spatial components of an ecosystem and their long-term interaction (Anderson and Battarbee 1994, Smol 2008). Lake sediments capture ecological conditions created by both gradual and rapid transitions in the environment and by both natural and anthropogenic stressors (e.g., Bracht-Flyr and Fritz 2012, Rühland et al. 2015). Compared to monitoring, paleolimnology faces its own challenges and limitations in that a record may integrate the impact of several interacting stressors, so teasing apart the relative importance of stressors over time requires multiple sedimentary proxies and other lines of evidence. There is potential for paleoecological research to address aspects of resilience theory, such as characterizing the shifting relative importance of environmental influences on an ecosystem over time.

This study provides a case study of the resilience of aquatic communities to fish stocking with specific consideration of the diatom (Bacillariophyceae) subfossil record. Carpenter et al. (1985) proposed the theory of trophic cascade interactions in food webs, in which the productivity of a lake ecosystem is regulated hierarchically, such that change at one trophic level due to biotic or abiotic factors effects changes in other trophic levels. Because lake productivity is at least partially regulated by abiotic factors, such as nutrient supply, the presence of a trophic cascade can account for differences in productivity among lakes with similar nutrient supplies yet different food web structure (Schindler 1978, Carpenter et al. 1985). A trophic cascade can be instigated by bottom-up or top-down influences. Bottom-up influences include abiotic conditions that affect primary productivity, such as nutrient availability and recycling, mixing, thermal stratification, and water clarity (Matson and Hunter 1992). Fish, as the apex predator in many freshwater ecosystems, exert a top-down influence in a trophic cascade through predation. Planktivorous fish often optimally forage on the largest available prey, resulting in a community dominated by smaller zooplankton, such as small cladocerans (Hall et al. 1976, Carpenter et al. 1985). Shifts in the zooplankton community, in turn, influence herbivory stress on phytoplankton, as well as nutrient cycling, which affects the community structure of the phytoplankton (Sommer 1989, Carpenter and Kitchell 1993, Attayde and Hansson 2001, Sarnelle and Knapp 2005). Thus, the presence or absence of fish reflects in the community composition of lower trophic levels, such as the phytoplankton (i.e., diatom) community.

Diatoms are reliable microfossil indicators of environmental characteristics (Smol 2008) and can respond to variation in fish population dynamics through fish impacts on nutrients and on food web structure (e.g., Finney et al. 2000). Alterations in trophic structure have been observed in lakes that have been purposefully stocked with fish for recreation or subsistence. In large lake systems, fish introductions can shift the food web to dominance by nonnative species, as occurred in Lake Superior (Kitchell et al. 2000); negatively stress native fish and macroinvertebrate species, as in Lake Victoria (Ligtvoet and Witte 1991); introduce preferential predation that can cause local extinctions and rearrange the food web, as in Flathead Lake (Spencer et al. 1999, Ellis et al. 2011); and alter the nutrient dynamics (Schindler et al. 2001). Trophic cascades have also been recorded in small shallow lakes, similar to those in this study. Fish introductions can alter the geographic distribution and biodiversity of invertebrates inhabiting naturally fishless lakes, such as those in Maine (Schilling et al. 2008, 2009) and the Canadian Rocky Mountains (Lamontagne and Schindler 1994). Previous studies have employed paleolimnological methods using diatoms to study trophic cascades. Linkages between diatom and cladoceran dynamics on grazing and nutrient recycling are evident in paleolimnological records (Harmsworth and Whiteside 1968, Velghe et al. 2012). Subfossil records of changes in the diatom assemblages in naturally fishless lakes were linked to the fish introductions through stocking programs in the Pacific Northwest in a manner similar to this study (Drake and Naiman 2000).

This study examines the effects of human-mediated fish stocking of lakes in Yellowstone National Park during the first half of the 20th century. Beginning as early as 1869, the U.S. Army stocked naturally fishless lakes with sports fish near their camps (Varley 1981). After the formal establishment of Yellowstone National Park in 1872, the National Park Service expanded the fish stocking program to include nearly all lakes in Yellowstone National Park. In 1955, the fish stocking program ended after certain species became invasive, and the service’s focus shifted to environmental conservation. Both the Army and the National Park Service kept detailed records of the species, number, and age of fish included in the stocking program (Table 1), allowing for a historical study of fish stocking impacts on naturally fishless lakes.

Table 1. Attributes of the lake study sites and their fish stocking histories.

	Lake of the Woods	Tanager Lake	Cascade Lake	Mallard Lake
Area (ha)	10.5	12.9	14.6	12.9
Maximum depth (m)	7	3.4	8.2	9.1
Transparency (m)	2.25	1.5	1.5	2.5
pH	6.44	6.69	7.7	6.76
TDS (g/L)	0.03	0.03	0.04	0.04
Conductivity (mS/cm)	0.04	0.04	0.06	0.03
Salinity (ppt)	0.02	0.02	0.03	0.01
Climate region	Summer-wet	Summer-dry	Summer-wet	Summer-dry
Hydrology	Snowmelt, runoff, no inlet/outlet	Seepage, spring, stream outlet	Springs, seepage, stream inlet and outlet	Snowmelt, seepage, no inlet/outlet
Bedrock	Surficial hotspring deposits, Obsidian Creek Member rhyolite	Surficial hotspring deposits, Lewis Canyon rhyolite	Hayden Valley rhyolite, Solfatara Plateau rhyolite, rhyodacite ashflow tuff, welded ash flows	Surficial detrital deposits, Obsidian Creek Member rhyolite-basalt mixed lavas
Catchment vegetation	Mossy bog, lodgepole pine forest	Boggy meadow, subalpine forest	Grassy meadow, pine forest	Lodgepole pine forest
Fish stock history	None	None	1928–1955	1919–1943
Number of fish per stocking event	n/a	n/a	50,000–550,000	4,000–120,000
Stocked fish species	n/a	n/a	Thymallus arcticus montanus, Oncorhynchus clarki	Oncorhynchus clarki, Salvelinus fontinalis
Current fish presence	Fishless	Fish	Fish	Fishless
Notes	Anoxic hypolimnion	Naturally has fish (minnows)

Note. Water chemistry parameters were measured on site while collecting sediment cores in 2013; these data are displayed along with climate region (Whitlock and Bartlein 1993), lake morphology and hydrology (Pierce 1987), bedrock (USGS 1972), and fish stocking records (Varley 1981). n/a, Not applicable.

Yellowstone National Park provides an intriguing site to study community resilience to the same event, fish stocking, because the lakes are spatially close yet are characterized by differences in microclimate. Specifically, Yellowstone is divided into 2 climate regions: The “summer-wet” region occurs in the northern portions of the park, where summer monsoonal rainfall from the Gulf of Mexico produces wet summers with comparatively dry winters (Whitlock and Bartlein 1993). In contrast, the “summer-dry” region occurs in the southwestern portions of Yellowstone, where the polar jet stream drives precipitation-bearing storms southward in winter, resulting in wet winters characterized by heavy snowpack that are coupled with relatively dry summers, because of subsidence associated with the Pacific subtropical high pressure system (Whitlock and Bartlein 1993). This spatial variation in climate provides the opportunity to examine how abiotic and biotic factors interact to influence an aquatic community’s resilience to a disturbance. This study considers the relative importance of trophic cascades in lakes from known fish stocking events in comparison to 20th-century climate trends. First, we hypothesize that if a trophic cascade did occur in a lake stocked with fish, the diatom assemblage would reflect the cascade effects via changes in its diatom assemblage and/or changes in algal productivity over time. Second, we hypothesize that the diatom assemblage of each lake will express resilience to a stressor if its assemblage remains constant over time, whereas a loss of resilience due to a stressor will result in a changing assemblage over time.

Methods

Study sites

All the study lakes are small, nongeothermal, subalpine, high-elevation lakes within the boundaries of Yellowstone National Park (Figure 1). We selected one pair of historically fishless lakes that were stocked with nonnative fish for sport. One lake (Mallard Lake) is currently fishless (stocked fish did not maintain viable populations), and one lake (Cascade Lake) is currently fish-maintaining (still inhabited by fish) at the time of the study. We also collected sediment cores from an additional pair of lakes (Feather Lake and Goose Lake) but these cores were not dated, so they are not discussed further in this article (see Supplemental Material). The dated replicate pair was used to assess the rates of response of the algal communities to fish introductions and the extent of trophic reorganization after the stocking program ended. We also selected 2 reference lakes that were never stocked with fish: Lake of the Woods is naturally fishless, and Tanager Lake is naturally inhabited by planktivorous minnows. In the subsequent text, the lakes are referred to by acronyms representing no stocking (NS) versus stocking (S), as well as the lake’s current status as either fishless (FL) or fish-maintaining (FM). The lakes were chosen for accessibility, similar environmental characteristics, and proximity to the paired lake in the replicate (Pierce 1987). The characteristics of each lake replicate pair are summarized in Table 1.

Figure 1. Sample site map. All lakes fall within the boundaries of Yellowstone National Park (thick red line). Lake paired replicates are designated by symbol (triangle, square, pentagon).

Mallard Lake was stocked with nonnative brook trout (Salvelinus fontinalis) and native cutthroat trout (Oncorhynchus clarki), mostly as either eyed eggs or yearlings (Varley 1981). Brook trout eat nearly any organism; they mostly predate upon aquatic insects, such as Ephemeroptera and Trichoptera. Larger adults will predate upon smaller fish, including their own young (Fry 1960, Brown 1971, Holton and Johnson 1996). Brook trout reach spawning age at about 2 years and have a life span of 5 years (Brown 1971, Holton and Johnson 1996). Native cutthroat trout are planktivorous as young and piscivorous as adults (Varley and Gresswell 1988, Gresswell 1995). Cutthroat trout reach spawning age at 5 years, and adults are larger (15 pounds) and longer lived (11 years) than brook trout (Varley and Gresswell 1988, Gresswell 1995). Cascade Lake was stocked with native Montana arctic grayling (Thymallus arcticus montanus), mostly as fry, advanced fry, or eyed eggs, except for one stocking event of blackspotted cutthroat trout (Oncorhynchus clarki) fingerlings in 1928 (Varley 1981). Grayling are generalist consumers of invertebrates; lake populations of grayling consume zooplankton as well as macroinvertebrates (Liknes and Gould 1987). Graylings have a fast growth rate and a life span of about 5 years (Liknes and Gould 1987).

Sediment sampling and analysis

A Pylonex HTH gravity corer collected a sediment core from the deepest part of each lake in July 2013. Core lengths were 22 cm for Cascade Lake, 26.25 cm for Mallard Lake, 34.75 cm for Lake of the Woods, and 36.5 cm for Tanager Lake. Cores were described and extruded in 0.25 cm increments in the field to preserve sediment for ²¹⁰Pb dating. The subsamples were individually packaged in Whirlpak bags and transported and stored in cold, dark conditions.

Approximately 0.5 g dry weight of sediment was subsampled from select core intervals of Lake of the Woods, Tanager Lake, Mallard Lake, and Cascade Lake. These sediment subsamples were dated using ²¹⁰Pb (Appleby and Oldfield 1978). The age model is based on a constant rate of supply (CRS) model, and the ²¹⁰Pb dates are converted to calendar ages using the Bayesian age model in the computer model Bacon (v. 2.2, Engstrom and Wright 1984, Blaauw and Christen 2011). Most sediment was young (<200 years) and dated within an error of 1–5 years.

Diatom preparation and analysis

For diatom analysis, approximately 0.25 g of sediment was first treated with 10% HCl to remove carbonates and then with 30% H₂O₂ to digest organic material. Samples were then rinsed 4 times with reverse-osmosis-purified water to remove supernatant fluids after treatment. Rinsed samples were dried onto coverslips and mounted on slides with Zrax, a permanent mounting medium with a high refractive index (Battarbee et al. 2001).

For community assemblage data, at least 300 diatom valves were counted on each slide (Battarbee et al. 2001). The assemblage data of dominant species (present in >5% abundance in at least one sediment sample) are shown in biostratigraphies as percent abundance relative to the total assemblage counted. The data are plotted against sediment depth as the primary axis and the ²¹⁰Pb age as a secondary axis. Zones defined by similar diatom assemblages are determined using constrained hierarchical clustering, in which clusters are constrained by sample order (i.e., depth interval). The clusters are determined by calculating the dissimilarity coefficient between data points using the CONISS incremental sum of squares method (Grimm 1987, Juggins 2015). Broken stick cluster analysis compares the hierarchical classification of the clusters and identifies the number of significant clusters, which have a cluster value higher than the broken stick value, as signifying temporal changes in the community structure that are identified graphically by horizontal zone boundaries (Birks and Gordon 1985, Bennett 1996). The hierarchical clustering and broken stick cluster analysis are calculated using the rioja package of R statistical software (v. 3.3.2; R Development Core Team 2014, Juggins 2015) and the biostratigraphies were plotted using the software package C2 (Juggins 2014). Although the sediment cores of Goose Lake (S/FM) and Feather Lake (S/FL) were not dated and hence are not discussed in the text, the diatom stratigraphies for these lakes are available in the Supplemental Material.

As a proxy of primary productivity, diatom concentration is calculated as a ratio to the presence of a known concentration of plastic microspheres (∼5 μm in diameter) and corrected to the dry weight of the sediment in the subsample. As a supporting proxy of primary productivity, stomatocysts of the algal group Chrysophyceae were enumerated, and their concentrations are calculated as a ratio to the known concentration of plastic microspheres and corrected to the dry weight of the sediment in the subsample. Similar to diatoms, chrysophytes are a flagellated, single-celled siliceous golden-brown algae that are often a dominant phytoplankton group in temperate, high-latitude, high-altitude lakes (Duff et al. 1995). The correlation between diatom and chrysophyte concentrations is calculated for each lake by linear regression to gauge how well chrysophytes serve as a supporting proxy for diatoms and general algal productivity.

A dissimilarity model based on squared chord distance (SCD) assessed shifts in the entire diatom assemblage during the sedimentary record as a proxy of resilience, similar to the method utilized by Drake and Naiman (2000), which also used diatom assemblages to evaluate restoration efforts in fishless lakes stocked with trout in Mt. Rainier National Park. SCD calculates a coefficient that dampens the effects of very rare taxa (noise), yet remains sensitive to the data signals of more common taxa (Overpeck et al. 1985). SCD can consider assemblages of varying diversities, because it is insensitive to species richness (Drake and Naiman 2000). The SCD (d_ij) is calculated for each depth interval using the equation Σ_k(P_ik^1/2 – P_jk^1/2)², in which p is the proportion of a given species (k) in the diatom assemblage. The baseline assemblage (j) is the proportion of all species present in the sum of samples prior to the beginning of fish stocking (prior to circa 1920), calculated specifically for each lake. The assemblage at a given time after circa 1920 (i) is calculated for each depth interval. Overpeck and others (1985) defined a benchmark of 0.15, meaning that if the SCD comparing 2 assemblages exceeds 0.15, then the 2 assemblages are considered to be different; this has been used to consider long-term assemblage changes in pollen (Overpeck et al. 1985) and diatoms (Flower et al. 1997, Drake and Naiman 2000). Using dissimilarity, we expect that a community assemblage that is highly resilient to a stressor would not become dissimilar enough to surpass an SCD of 0.15.

Dissimilarity trends are evaluated relative to patterns predicted by the 3 hypotheses (Figure 2(a, b)). First, a dissimilarity time series that shows stability around the pre-fish stocking baseline and does not vary above the critical value of 0.15 suggests high resilience to the stress of fish stocking. Second, a dissimilarity time series that shows a loss of stability after the introduction of fish (i.e., dissimilarity values increase above the critical value of 0.15) but regains a baseline suggests a trophic cascade effect temporally associated with fish stocking. Within this, the assemblage is expected to either recover the original baseline value or establish a new baseline value (i.e., alternate stable state). Third, a continually increasing dissimilarity over time suggests a loss of resilience and potential movement toward a regime shift, particularly if the trend does not plateau. If the trend is not temporally associated with either the initiation or conclusion of fish stocking, the trend is likely being driven by a stressor other than fish stocking.

Figure 2. Dissimilarity analysis. (A) Hypothetical trends of diatom assemblage response to trophic cascade effects in terms of resilience, displayed as (B) corresponding dissimilarity trends that would reflect resilience or regime shift as changes in the community assemblage. The calculated dissimilarity trends are shown for the entire diatom assemblage of (C) Lake of the Woods, (D) Tanager Lake, (E) Cascade Lake, and (F) Mallard Lake. The red dotted line represents a critical value of 0.15; points that fall above this line have significant dissimilarity between the baseline assemblage (1850–1920) and the assemblage of the plotted time interval.

Results

Sediment dating

The 4 dated lakes showed variation in their sediment ²¹⁰Pb activity profiles (Table 2, Supplemental Material), ²¹⁰Pb-based dates and dating errors (Figure 7), and sediment accumulation rates (Table 2, Supplemental Material). Lake of the Woods has a near-exponential ²¹⁰Pb activity profile, with the oldest reliable date of 1812 ± 22 years and a low sediment accumulation rate of about 0.09 cm/yr since circa 1900, with a broad plateau of sediment accumulation from circa 1860 to 1940. Tanager Lake has an irregular ²¹⁰Pb activity profile, with the oldest reliable date of 1814 ± 27 years. Tanager Lake had a moderate sediment accumulation rate of about 0.17 cm/yr since circa 1900, with 2 peaks in accumulation circa 1920–1950, and rising accumulation rates after circa 1990. Mallard Lake also has an irregular ²¹⁰Pb activity profile, with the oldest reliable date of 1834 ± 12 years, a moderate rate of sediment accumulation of about 0.13 cm/yr since circa 1900, a peak accumulation rate from circa 1910–1950, and a suspected shoreline slump circa 1950. Cascade Lake has a monotonic down-core decrease of ²¹⁰Pb activity with an oldest reliable date of 1856 ± 16 years. Cascade Lake had a moderate sediment accumulation rate of about 0.13 cm/yr since circa 1900, with a small peak in accumulation circa 1920, and rising accumulation rates after circa 1990.

Figure 3. Diatom biostratigraphy of Lake of the Woods, Yellowstone National Park, WY, since 1850. The dotted line designates a zone of significant change in the diatom assemblage as identified by cluster analysis shown in the dendrogram. Diatom and chrysophyte concentrations provide proxies of primary productivity. Lake of the Woods was never stocked and is naturally fishless.

Table 2. ²¹⁰Pb parameters and sediment accumulation rates for cores from Yellowstone National Park, including 2 reference lakes and 2 lakes included in the fish stocking program.

Coring Site	Lake of the Woods	Tanager Lake	Mallard Lake	Cascade Lake
Cumulative unsupported ^210Pb (pCi/cm^2)	13.2	39.1	25.8	28.4
Unsupported ^210Pb concentration at surface (pCi/g)	61.8	37.2	77.9	29.8
Supported ^210 Pb (pCi/g)	1.03 ± 0.03	1.49 ± 0.09	0.97 ± 0.03	0.86 ± 0.06
Number of supported samples	4	6	6	4
Mean sedimentation rate since 1900	0.09	0.17	0.13	0.13
Mean ^210Pb flux (pCi/cm^2-yr)	0.42	1.31	0.83	0.95

Trends in diatom assemblage

All sediment cores displayed a viable diatom subfossil record. We present the results by paired replicate. Although the records dated to ages varying between 1812 and 1856, only the results for the mutual time period present in all records (since 1850) are presented, to aid in clear discussion. This provides a baseline of ∼30 years before the fish stocking program in Yellowstone commenced.

Lake of the Woods and Tanager Lake (never stocked)

Lake of the Woods (NS/FL) is a naturally fishless lake that was never stocked. The diatom assemblage had a high diversity of the benthic raphid genus Pinnularia, although most species were rare. The record over ∼163 years underwent a general shift in dominance from the planktic species, Aulacoseira alpigena-distans (Grunow) Krammer 1991, to benthic species (Figure 3). The shift began at circa 1935, with a peak in abundance of 2 species of the benthic genus Pinnularia, followed by several other benthic or tychoplanktic species (attached taxa that can survive when entrained into the plankton). There is a general trend of decreasing diatom concentration, which we use as a proxy for primary production by this algal group. The chrysophyte statocyst concentration, which we use as a supporting proxy for primary production of other algae, remains relatively constant through time. Diatom and chrysophyte concentrations were not correlated (r² = 0.03, P > 0.05, α = 0.05); a notable exception is a brief concurrent peak in both diatom and chrysophyte concentration at circa 1972.

Figure 4. Diatom biostratigraphy of Tanager Lake, Yellowstone National Park, WY, since 1850. The dotted lines designate zones of significant change in the diatom assemblage as identified by cluster analysis, shown in the dendrogram. Diatom and chrysophyte concentration provide proxies of primary productivity. Tanager Lake naturally has fish and was never stocked.

Tanager Lake (NS/FM) was never stocked but has a natural population of small planktivorous minnows. Tychoplanktic and benthic species predominate in the assemblage (Figure 4). No single species strongly dominates the fossil diatom assemblage at a given time. Notable assemblage shifts occurred at circa 1878, with a decrease in Diatoma moniliformis Kützing 1833, and in circa 1956, when Pseudostaurosira trainorii Morales 2001, Staurosira construens Ehrenberg 1843, and Staurosirella pinnata (Ehrenberg) Williams and Round 1987 increased. The planktic diatom Aulacoseira granulata (Ehrenberg) Simonsen 1979 also increased around this time. Diatom and chrysophyte concentrations were correlated (r² = 0.82, P < 0.05, α = 0.05).

Figure 5. Diatom biostratigraphy of Mallard Lake, Yellowstone National Park, WY, since 1850. The dotted lines designate zones of significant change in the diatom assemblage as identified by cluster analysis, shown in the dendrogram. Diatom and chrysophyte concentration provide proxies of primary productivity. Fish stocking numbers reflect National Park Service records of stocking with native cutthroat trout (Oncorhynchus clarki) and nonnative brook trout (Salvelinus fontinalis); the background shading marks the stocking program duration at this lake (Varley 1981). Originally fishless, Mallard Lake returned to a fishless state after being stocked with fish for 24 years.

Mallard Lake vs. Cascade Lake (stocked)

Mallard Lake (S/FL) was stocked with fish between 1919 and 1943 and has since returned to a fishless state. Over ∼163 years, the diatom assemblage showed a general trend of decreasing planktic species abundance and increasing benthic species abundance (Figure 5). Aulacoseira alpigena-distans dominated until circa 1872; concurrently, planktic Fragilaria tenera (W. Smith) Lange-Bertalot 1980 and tychoplanktic S. pinnata and Pseudostaurosira brevistriata (Grunow) Williams and Round 1987 decreased. Between circa 1872 and circa 1962, A. alpigena-distans decreased and was replaced by Aulacoseira ambigua (Grunow) Simonsen 1979 and the benthic species Navicula cryptocephala Kützing 1844, Sellaphora rexii Potapova and Ponader 2008, and Stauroneis pikuni Bahls 2010. After circa 1962, tychoplanktic species returned to dominance after almost a century of low abundance (<4%), in particular P. brevistriata and S. pinnata. Diatom concentration was low until circa 1979, when it increased. Chrysophyte concentration has a similar increasing trend over the past 40 years (Figure 5) and is correlated with diatom concentration (r² = 0.76, P < 0.05, α = 0.05). Neither diatom nor chrysophyte productivity increased during the stocking period.

Figure 6. Diatom biostratigraphy of Cascade Lake, Yellowstone National Park, WY, since 1850. The dotted lines designate zones of significant change in the diatom assemblage as identified by cluster analysis, shown in the dendrogram. Diatom and chrysophyte concentration provide proxies of primary productivity. Fish stocking numbers reflect National Park Service records of being stocked with native Montana arctic grayling (Thymallus arcticus montanus) and native cutthroat trout (Oncorhynchus clarki); the background shading marks the stocking program duration at this lake (Varley 1981). Originally fishless, Cascade Lake independently continues to maintain a fish population after being stocked with fish for 27 years.

Cascade Lake (S/FM) was stocked with fish between 1919 and 1943 and continues to maintain a fish population after stocking ended. Over the past ∼150 years, the record showed a general trend of decreasing planktic species abundance and increasing benthic species abundance (Figure 6). In circa 1944, the dominant species A. alpigena-distans, A. ambigua, A. granulata, Lindavia bodanica (Eulenstein ex Grunow) Nakov et al. 2015, and Discostella pseudostelligera (Hustedt) Houck and Klee 2004 decreased and were replaced by Asterionella formosa Hassall 1850, Tabellaria fasciculata (Agardh) Williams and Round 1986, and P. brevistriata. Another shift at circa 1983 resulted in dominance by P. brevistriata, F. tenera, S. pinnata, and Fragilaria crotonensis Kitton 1869. Diatom and chrysophyte concentrations are correlated (r² = 0.89, P < 0.05, α = 0.05); both peaked briefly at circa 1960, then increased after circa 1983.

Figure 7. Sediment age (year C.E.) based on sedimentary ²¹⁰Pb activity (pCi/g) of Lake of the Woods, Tanager Lake, Mallard Lake, and Cascade Lake.

Trends of assemblage change

All 4 lakes show a trend of increasing dissimilarity in the diatom assemblage from 1850 to 2013. Lake of the Woods (NS/FL), Tanager Lake (NS/FM), and Cascade Lake (S/FM) show a trend that increases above the critical value of 0.15 beginning at circa 1920, which is within a decade of the fish stocking program, though only one of the lakes was stocked. Mallard Lake (S/FL) shows the trend beginning earlier than the stocking program in circa 1910. Mallard Lake (S/FL) shows a new, although highly variable, baseline (i.e., alternate stable state) during the period of fish stocking from circa 1925 to 1955, then an increasing dissimilarity during the recovery period to a fishless state.

Discussion

The trend of increasing dissimilarity becomes pronounced in all lakes at circa 1955, within a decade of the end of the fish stocking program. Mallard Lake does not show a recovery to the original fishless community assemblage, and the pattern suggests a continual loss of resilience that has not yet stabilized. Lake of the Woods (NS/FL) and Cascade Lake (S/FM) show similar trends of increasing dissimilarity. Tanager Lake (NS/FM) shows a muted trend compared with the other 3 lakes; its diatom assemblage has stayed more similar between the past and present than the other 3 lakes. The similar trends in the lakes despite different stocking histories and different current fish status, combined with the lack of temporal coordination of the trend with fish stocking (with the possible exception of Mallard Lake), suggest that an environmental influence other than fish stocking is driving the loss of resilience.

The basin morphology of some lakes was insufficient to support fish populations. Adult brook trout and cutthroat trout tend to spawn in gravel-bottomed small streams, where they construct redds (Brown 1971, Varley and Gresswell 1988, Holton and Johnson 1996). Mallard Lake (S/FL) has no stream inlet or outlet for fish to spawn, which is likely why it failed to sustain a fish population after the stocking program ended. In contrast, Cascade Lake (S/FM) has a stream inlet and outlet, providing the necessary spawning habitat required to maintain a fish population after the stocking program ceased. Moreover, trout fry seek shallow waters and littoral vegetation for shelter, as well as macroinvertebrate prey while they develop (Brown 1971, Varley and Gresswell 1988, Holton and Johnson 1996). Mallard Lake (S/FL) is situated in lodgepole pine forest on a volcanic dome, resulting in a basin with steep or slumped shorelines and limited littoral habitat. In contrast, Cascade Lake (S/FM) is located in grassy meadow terrain with vegetation extending up to the shoreline of the lake, producing more suitable nursery habitat for fry.

In comparing the paired replicate of stocked lakes, Mallard Lake (S/FL) displayed a trend of trophic cascade impact most similar to that hypothesized, in that the dissimilarity of diatom assemblages increased, and the lake established a new baseline during the stocking period (Figure 2). This baseline includes high variability, which may be due to the limited but frequent effects of large fish stocking events in which large numbers of fish were added at intervals but could not successfully reproduce, potentially resulting in a fish–fishless–fish–fishless cycle every few years. This new baseline ends at circa 1955, at 12 years after the final stocking event in Mallard Lake and perhaps when Mallard returned to a fishless state. From 1955 onward, the dissimilarity trend shows a loss of resilience likely not associated with fish. In contrast, Cascade Lake (S/FM) displays small increases in dissimilarity during fish stocking compared to the baseline (Figure 2), suggesting that the fish, even though the population was maintained and had the capacity to produce larger trophic effects, did not exert a large influence on the diatom community.

The observed shift of diatom assemblages in both Mallard (S/FL) and Cascade (S/FM) to benthic species dominance could be due to trophic cascade effects. The stocked fish species are known to feed upon macroinvertebrates and zooplankton, and trout in particular will feed upon benthic macroinvertebrates and, by so doing, alter nutrient cycling (Hecky and Hesslein 1995) and sever trophic connections in watersheds reliant upon larval and adult aquatic insects (Finlay and Vredenburg 2007). The planktivorous diet by stocked fry would relieve herbivory pressure on benthic diatoms, resulting in an increased population; this type of top-down control has been shown to be important in shallow lakes due to the high biomass of benthic invertebrates and fish relative to depth (Jeppesen et al. 1997). Diatoms and chrysophytes, as proxies of primary productivity, did not indicate a significant impact of trophic structure (fish stocking) on primary production. The rhyolitic bedrock of much of Yellowstone provides high external loading of P to these regional lakes (Theriot et al. 1997), perhaps alleviating P limitation in most systems, such that any altered cycling from a trophic cascade would have little effect on the diatom community. Benthic species obtain a substantive portion of their nutrients from the sediment rather than nutrient cycling in the water column (Engstrom et al. 1991), so the high abundance (>80%) of benthic species in all of the lakes may mute any signal of trophic cascades on nutrient cycling.

The record of trophic cascade effects in Yellowstone lakes is muted or not evident, and stocked lakes show no evidence of recovery after stocking ended. Similarly, the diatom communities of Mt. Rainier National Park demonstrated a response to trout introduction, yet also failed to recover to baseline conditions as a fishless community after the cessation of fish stocking on time scales similar to the Yellowstone lakes in this study (Drake and Naiman 2000). Overall, the timing of the shifts and the shared trends among lakes do not support the hypothesis of trophic cascade as the most parsimonious explanation for changes in the diatom assemblages. Instead of trophic cascade, the weight of evidence suggests that the lakes recorded diatom assemblage shifts in response to some other factor, such as changing climate. In addition to the shift to benthic diatom assemblages in Mallard (S/FL) and Cascade (S/FM), a similar trend was observed in Lake of the Woods (NS/FL), which one would not expect to observe if the trend is explained solely by trophic cascade effects (Figure 8). Moreover, all 4 lakes, regardless of stocking history or current fish population status, displayed similar trends of increasing dissimilarity over the past 70–100 years (Figure 2), suggesting a widespread environmental influence. Shifts to benthic species dominance may be interpreted as a response to lower lake level, increased water clarity, and/or nutrient limitation, such as P or Si (Engstrom et al. 1991).

Figure 8. Climate parameters and diatom benthic community assemblage in Yellowstone National Park since 1850. Minimum temperature (a proxy of winter temperature), maximum temperature (a proxy of summer temperature), average annual temperature (C), and average precipitation (mm) with a 5-year smoothing spline for the past ∼100 years from the Mammoth Hot Springs weather station. Snow water equivalent (a proxy for snowpack) with a 20-year cubic smoothing spline for the past ∼200 years is shown relative to the observed 1 Apr snow water equivalent (cm) observed in northern Yellowstone (Pederson et al. 2013). Trends in benthic diatom assemblage dominance are shown in terms of percent relative abundance of all benthic species present in the unique assemblages of 4 lakes.

Lower lake level increases the amount of littoral habitat and hence the exposure to light for benthic diatoms or the macrophytes on which they live (Engstrom et al. 1991). The increases of benthic or epiphytic fragilarioid species in all of the lakes in this study, including S. pinnata, S. construens, and P. brevistriata, have been noted in other small lakes in the northern range of Yellowstone from 1920 to 1960, in probable response to lower lake level or oligotrophic conditions (Engstrom et al. 1991). Lower lake level may have resulted from declining snowpack in the region. Pederson and others (2013) modeled long-term temperature and precipitation data from the region to assess periods of high or low snowpack, which strongly influences regional water availability in the following summer (Figure 8). The northern Rockies had high snowpack in the periods of 1900–1928, 1942–1957, and 1964–1980, and low snowpack in the periods of 1929–1941 and 1958–1963. Since 1980, the region has had persistent declines in snowpack, coupled with warmer earlier springs.

In Lake of the Woods, the shift to benthic diatom dominance began at circa 1935 during a period of low snowpack. The trend toward benthic diatom dominance in Cascade Lake (circa 1945) and in Mallard Lake (circa 1965) began during periods of high snowpack that immediately followed periods of low snowpack, and the trend intensified in both lakes after circa 1980, the beginning of persistent snowpack decline. Lake of the Woods and Cascade Lake are located in the region of Yellowstone that receives the majority of its precipitation in summer months (the “summer-wet” region). In this region, lake level not only is influenced by summer convective storms but also is affected by winter snowpack and the associated meltwater in spring to early summer (Whitlock et al. 2012). Although it is in the “summer-dry” region, Mallard Lake displayed a trend similar to the lakes in the summer-wet region, likely because it is also fed by snowmelt and seepage with no inlet or outlet (Pierce 1987). Declining snowpack results in less runoff and less effective moisture in the spring and early summer, leading to lower lake levels. The lowered runoff also reduces the external load of nutrients and organic material to lakes during the spring snowmelt, to the deficit of the growing season’s phytoplankton.

In contrast, Tanager Lake and, as a note, Goose and Feather lakes (see Supplemental Material) did not have a trend toward benthic dominance but instead were dominated by benthic species throughout the record. These lakes are located in the “summer-dry” climate region of Yellowstone, which does not receive substantive precipitation in the summer, and are fed by groundwater and springs. Thus, the summer-dry lakes may be more resilient to declines in snowpack if they are fed by groundwater. Instead of lake level, summer-dry lakes may be dominated by benthic species as a result of other factors. Tanager Lake has an extensive shallow littoral area that provides excellent habitat for benthic growth. It is also located in a marshy meadow instead of tall forest, so it is possible that the lake is subject to high wind mixing, which favors tychoplanktic species, such as those dominant in Tanager Lake.

Overall, the diatom assemblage changes are consistent with independent evidence of recent regional climate change toward shorter winters and the consequent impact on lake physical structure, including a reduction in the length of ice cover, earlier and prolonged spring water-column mixing, longer summers with stable thermal stratification, and changes in water level and water clarity caused by reduced snowpack. This suggests the importance of physical factors on biological structure via their influence on niche habitats and nutrient cycling (Figure 9, Rühland et al. 2008, 2015, Fritz and Anderson 2013).

Figure 9. A generalized conceptual flow chart of the hierarchy of environmental controls on the aquatic ecosystems, modified from Fritz and Anderson (2013) to represent variables considered in this study. The hierarchy is represented by the gradation of background color. The thick arrows represent interactions among environmental drivers. Thin arrows represent interconnections between the 2 main categories of environmental controls, climate, and catchment. The fish stocking/trophic cascade entry is characterized by a dotted line to denote its hypothesized rank inthe hierarchy in this study. Note that this flow chart only depicts environmental variables considered in this study and is not comprehensive.

This paleolimnological study supports the hypothesis that climate often provides a first-order control on the physical and biological structure of small, high-elevation, temperate lakes with little direct human impact (Fritz and Anderson 2013). This influence appears to function on varied time scales, from decades to centuries, and at various spatial scales, from the microclimate of a watershed to regional or hemispheric climatic trends (Figure 9). Results from this case study suggest climate influenced diatom assemblages through less precipitation, lower lake levels, and changes in lake thermal properties such as spring mixing and summer stratification. The climatic mechanisms seem to either override or compound the effects of trophic cascades on diatom assemblages, because as phytoplankton, diatoms respond along a species-specific sensitivity gradient to various environmental stresses that influence both assemblage and productivity (Smol 2008). This highlights an important caveat to consider in paleolimnological studies with a potential application to lake management: Factors can occur at different scales at the same time, potentially resulting in similar signals. For this reason, it is important to consider paleolimnological data holistically using multiple lines of evidence.

Biological interactions, such as predation, herbivory, and productivity, were not major influences on the lakes’ biological structure over the time scales considered in this study. This is surprising, considering evidence from other Northern Hemisphere lakes showing clear changes in response to trophic cascades (e.g., Jeppesen et al. 2001). The aquatic communities may have been buffered by aspects of their environment, such as nutrient-rich bedrock that relieves dependence on nutrient recycling, or by the failure of many fish species to establish permanent populations. It is likely that in this particular case, the phytoplankton communities responded more strongly to the bottom-up influences of thermal structure and lake level than to top-down controls from trophic cascades. This interpretation is corroborated by evidence from Yellowstone Lake, a much larger lake in the region, where the introduction of lake trout exerted a large negative influence on the native fish population and lower trophic levels (Dunham et al. 2004, Tronstad et al. 2010), yet the paleolimnological record shows a relatively stable diatom assemblage through time (Kilham et al. 1996, Interlandi et al. 1999).

Relevance to lake management

Resilience is an increasingly important consideration in the development of management protocols (e.g., Folke 2003), yet resilience is a complex concept that is still largely abstract. Recent attempts to create a quantitative framework to assess resilience are promising (e.g., Angeler and Allen 2016, Baho et al. 2017), but the exact methodologies are still under development. In freshwater systems, paleolimnological data offer the benefit of a long temporal scale with which to identify thresholds and regime shifts (e.g., Spanbauer et al. 2014). However, although diatoms are effective as a proxy of environmental drivers, such as climate or system productivity (e.g., Lu et al. 2017), they may provide an incomplete picture of smaller-scale dynamics within the system, including ecosystem responses to trophic cascades. In these situations, additional proxies, such as zooplankton (e.g., Jeppesen et al. 2001, Knapp et al. 2001) and chironomids (e.g., Lamontagne and Schindler 1994), may provide a more complete picture of the entire system. Thus, lake management frameworks seeking to assess the resilience of the system using paleolimnological data may require different proxies or a combination of proxies depending on the environmental driver or the scale of change in question. In this case study of Yellowstone, the diatom assemblages displayed high resilience to the localized and short-term stress of fish introductions, but displayed a lower resilience to the larger spatial and longer temporal scale stress of climate change.

Resilience is, by Holling’s (1973) definition, the ability to not change critical structure or function in response to stress. From this perspective, aquatic communities in these small, high-elevation lakes were resilient to fish stocking in all cases, but not resilient to climate change in some cases. This type of knowledge can inform management by guiding consideration of the type of resilience that is the goal of a given management plan. For example, the management plan with an underlying goal of functional resilience, which can be supported by species redundancy as is likely the case in the Yellowstone diatom assemblage, would likely look different from a plan to promote resilience through high biodiversity of rare or highly specialized species (Baho et al. 2017). Moreover, it is important to frame any management plan in the relative terms of “resilience of what to what” (Carpenter and Kitchell 1993). If a community displays resilience in one aspect of the system, an adaptive management plan can focus its effort and finances on bolstering a less resilient aspect of the system in order to achieve greater resilience of the overall system.

Panarchy is closely related to the concept of a hierarchy of environmental controls in that changes are often effected by higher order mechanisms. For example, climate influenced changes in lake level, which in turn resulted in changes in the diatom assemblage to a dominance of benthic taxa. Even so, panarchy suggests that lower order controls may accumulate to accomplish higher order change or to mitigate higher order influences. For example, it is possible that nutrient-rich bedrock and lake morphology mitigated the impacts of a higher order control, such as nutrient cycling and trophic structure. This presents 2 important considerations for lake management. First, lower order environmental controls may bolster resilience in a community to the effects of higher order environmental stress. In this case study, water level as determined by hydrological budget seems important to the resilience of small, shallow, high-elevation lakes to a high-order control of climate change. Second, the classification of lower orders of control in a hierarchy may be unique to each system of study. Thus, understanding of the hierarchy of environmental controls for an ecosystem, and the interplay among controls, may need to be site specific and continually adapted through time in order for effective adaptive management applications.

The 4 lakes in this study showed varying trends resulting from their unique diatom assemblages and locations. A general increase in fragilarioid species and an increase in benthic diatom species overall suggest lower lake levels, possibly in relation to declining snowpack and warmer temperatures. Other trends in specific diatom species also may result from changes in ice cover, seasonality, thermal structure, and wind mixing of the water column. Thus, in oligotrophic, relatively pristine high-elevation lakes in the Yellowstone region, trophic cascade impacts from fish stocking appear to have had negligible influence on the phytoplankton community in comparison to other anthropogenic influences, such as climate warming. In Yellowstone, the diatom communities show evidence of climate-driven shifts in community structure since circa 1910 to 1920, increasing after circa 1955. As of 2013, the communities have not reorganized in an alternate stable state, suggesting continued change for the foreseeable future. In terms of the management of such small, shallow, high-elevation lakes, these patterns suggest that adaptive management plans need to focus on actively adaptive systems.

Acknowledgments

Dr. Robert Gresswell and Dr. Cathy Whitlock provided valuable insight into study site selection. We also thank Dr. Gresswell for providing historical records of fish stocking and giving valuable feedback on the article. Dr. D. Marie Weide, James Benes, and Dave Mosicki assisted with field sampling. Sediment ²¹⁰Pb dating was done under the supervision of Dr. Daniel Engstrom at the St. Croix Watershed Research Station. Snow water equivalent data are courtesy of Dr. Greg Pederson. The site map was produced by Dr. Xavier Benito.

Additional information

Funding

This research was funded by National Science Foundation 10.13039/10000 (NSF) grant EAR-0816576 to S. Fritz and an NSF-IGERT fellowship (0903469), Geological Society of America10.13039/100005720 grants, and a Paleontological Society10.13039/100010841 grant to V. Chraïbi.