A Holocene landscape dynamic multiproxy reconstruction: How do interactions between fire and insect outbreaks shape an ecosystem over long time scales?

At a multi-millennial scale, various disturbances shape boreal forest stand mosaics and the distribution of species. Despite the importance of such disturbances, there is a lack of studies focused on the long-term dynamics of spruce budworm (Choristoneura fumiferana (Clem.)) (SBW) outbreaks and the interaction of insect outbreaks and fire. Here, we combine macrocharcoal and plant macrofossils with a new proxy—lepidopteran scales—to describe the Holocene ecology around a boreal lake. Lepidopteran scales turned out to be a more robust proxy of insect outbreaks than the traditional proxies such as cephalic head capsules and feces. We identified 87 significant peaks in scale abundance over the last 10 000 years. These results indicate that SBW outbreaks were more frequent over the Holocene than suggested by previous studies. Charcoal accumulation rates match the established fire history in eastern Canada: a more fire-prone early and late Holocene and reduced fire frequency during the mid-Holocene. Although on occasion, both fire and insect outbreaks were coeval, our results show a generally inverse relationship between fire frequency and insect outbreaks over the Holocene.


Introduction
Boreal forests are subjected to various natural disturbances that operate at different spatiotemporal scales [1,2] affecting biodiversity and ecosystem dynamics [3,4]. In Canada, fire and insect outbreaks are the two main disturbances that shape forest stand age and composition [5,6]. Although anthropic disturbances (forest harvesting, hydroelectrical infrastructures, oil and gas extraction, etc.) have had a significant impact on terrestrial ecosystems over the last centuries [7], natural disturbances continue to operate at a pluri-millennial scale [8,9].
Dendrochronological studies, through use of fire scars and reductions in tree ring growth, provide valuable insight into the history of fire and insect outbreaks. However, although forest compositions in the context of Holocene climate variation. We use lepidopteran wing scales in combination with charcoal and plant macrofossils (e.g. needles, leaves, seeds) to reconstruct, for the first time, a Holocene-scale history of the boreal landscape that includes changes in vegetation, fire frequency, and insect outbreaks.
We hypothesize that lepidopteran scales should be found throughout the entire stratigraphy. Scattered episodes of high scale abundance should correspond to periods of warmer temperatures and lower fire activity which contribute to a good development of the insect populations and a high abundance of mature host trees. Scale abundance should also be correlated with the presence of fir and spruce macrofossils within the sediment record. In a much shorter time scale, higher scale abundance should also increase short-term fire frequency as insect-related defoliation increase the risk of ignition, providing an important quantity of flammable fuel to the forest ground [25].

Materials and methods
The study was conducted in the university experimental forest and we did not need any permission to sample the sites for this study. The field study did not involve endangered or protected species. This study did not imply vertebrate species.

Site selection and sediment coring
The study area, Lake Flévy (48˚13'00"N; 71˚12'57"W), is in southeastern Quebec at the interface between the eastern balsam fir-white birch domain and the western balsam fir-white birch domain (Fig 1). The former domain is known to be drier and to has a shorter fire return interval [26]. Lake Flévy is also located in the northern portion of the present-day spruce budworm distribution providing possible insight on the evolution of species distribution under changing climate conditions [27]. The surrounding landscape is characterized by highly sloping, forestcovered hills dominated by spruce-fir stands, particularly on well-drained mesic to hydric locations (hill tops and valley bottoms). Lake Flévy is surrounded by even-aged stands of trembling aspen (Populus tremuloides (Michx.)) and mixed even-aged stands of black spruce and aspen. These forest stands originated after an intense fire that occurred in 1922 AD [28] (Fig 1).
SBW has been present in the area for a minimum of 8 240 years [9] and a number of outbreaks have occurred over the last three centuries, becoming more frequent during the 20th century [29]. Aerial surveys indicate that since 2012 AD, SBW-related defoliation has occurred in the resinous stands around the lake reaching severe levels after 2015 AD [30,31] (Fig 1). Regional climate is continental subpolar with a mean annual temperature of 0˚C, mean annual precipitation of 1000 mm and an average growing season of 155 days [32].
The study lake is located in the university experimental forest of Simoncouche and no permission was required in order to sample this site. The field study did not involve endangered or protected species.
The lake was chosen for its small size (2.33 ha) and limited inflow and outflow to the lake, thereby favouring a relatively high sedimentation rate and the retention of a local signal within the sediment record (Fig 1). Two cores were extracted from the deepest portion of the lake (4.8 m) using a Livingstone piston corer. The first core (AG1) was analyzed at the Université du Québec à Montréal (UQAM) for plant macrofossils and charcoal. The second core (AG2) was analyzed at Université du Québec à Chicoutimi (UQAC) for SBW microfossils.

Sample preparation
The chronological framework for both cores was determined using AMS radiocarbon dating of organic sediments. Samples from AG1 were analyzed at the Beta Analytic Lab in Miami, Florida. Samples from AG2 were prepared in the radiochronology laboratory of Université Laval's Centre for Northern Studies, then sent to the Keck Carbon Cycle AMS Laboratory, University of California, Irvine. The InCal13 database was used for calibrating 14 C dates [33]. Sediment accumulation rates were calculated applying a 3 rd degree polynomial model using the Clam 2.2 software (Fig 2). All dates are expressed as calibrated years BP (cal yr BP) ( Table 1).
Both cores were subsampled at a centimeter-scale resolution. For charcoal extraction, a ca. 1 cm 3 subsample was collected at each centimeter of AG1. This sediment was deflocculated in a 100 mL, 3% NaP 2 O 4 solution for 3 h. Each sample was then sieved through a 160 μm mesh to collect charcoal derived from local fires [34], then soaked in a 10% NaOCl solution to bleach the organic matter and help discriminate charcoal fragments.
The remaining AG1 sediment was deflocculated and then sieved through a 160 μm mesh. The recovered >160 μm fraction was transferred to petri dishes. Macrofossil remains, such as needles, twig, seeds, and cone scales, were analyzed using a binocular microscope and identified based on comparisons with reference collections [35][36][37].
Samples from AG2 were prepared following Navarro et al. [24]. Each sample was dried at 105˚C for 24 h to recover a 0.5 g subsample of dry sediment (± 5 cm 3 ). These samples were then heated in a 100 mL 10% potassium hydroxide (KOH) solution at 70˚C for 30 min or until complete deflocculation had occurred. The slurry was then sieved through a 53 μm mesh to retain most scales. We centrifuged the samples at 500 rcf for 10 min in a 10 mL sucrose solution (relative density = 1.24) to remove higher density particles. The centrifuging was repeated three times. After each run, we recovered the supernatant, refilled the vial with the sucrose solution, and centrifuged again. We combined the three supernatants in a 50 mL plastic vial and, to precipitate scales and any remaining particles, centrifuged the combined supernatant at 3900 rcf for 20 min. The final pellet was mounted onto microscope slides for microfossil counting.

Data handling and identification of peaks
Accumulation rates were calculated using CharAnalysis 1.1 [38]. In order to balance sedimentation rates variations amoung the two cores, each proxy concentration values and deposition times were interpolated to pseudo-annual intervals. The resulting values were then integrated over 5-year intervals and divided by the average deposition time over those intervals [39]. The influx series so obtained (C int ) allows for core to core comparaision and peak analysis. A stratigraphy was developed using the accumulation rates of each indicator. In Psimpoll 4.27, we used stratigraphically constrained cluster analysis (CONISS) to define three different assemblage zones [40]. Principal component analysis (PCA) was used to identify the main variables that explain the evolution of the ecosystem within these zones.
To extract fire events from the charcoal stratigraphy, a background component (C back ) was defined using a Lowess smoothing that was robust to outliers and that had a smoothing window of 500 years. C back have been extracted from the interpolated serie of raw data (C int ) to define a peak serie (C peak ) as a residual of C int -C back . A noise component (C noise ) was calculated for the C peak series using a Gaussian mixture model. A threshold for fire detection (C fire ) was defined based on the 99 th percentile of the C noise distribution. Cfire samples and samples preceeding the Cfire sample that have a >5% probability of being from the same Poisson distribution were discarded. From the fire event dates, we calculated fire frequency over a 500 years bandwidth and smoothed this using a Lowess smoother.
This procedure was also used for insect outbreak detection and frequency. In the same way that a charcoal peak is associated with a fire event, a scale peak is associated by charanalysis with a period of high lepidotera abundance. It is important to note, however, that at this temporal scale, a peak could represent multiple outbreaks. All scales from samples that were associated with a peak were identified so as to confirm a correspondence to a known outbreaking lepidopteran morphotype [24]. ANOVA was used to determine the mean differences between the relative occurrence of each species morphotype, and a Tukey's range test determined whether the means were significantly different.

Results
Core AG1 was 519 cm long and was dated to more than 6000 cal yr BP (Fig 2). One date was excluded from the age-depth model as it was considered as outlier (Table 1). Indeed, the radiocarbon date at the bottom of AG1 (Beta-401641) was probably influenced by the presence of stromatolites in the sampling lake (B. Lapointe, personal communication), probably originating from the lac Mistassinni [41] The sediment accumulation rate was relatively stable across the core varying from 6 to 16 yrÁcm -1 . Core AG2 was 511 cm long and represented almost 10 000 yr of deposition. The sediment accumulation rate of AG2 was more than two times slower than that of AG1 from 10000 to 6000 cal yr BP, then was more similar to that of core AG1 over the last 6000 years (Fig 2).
The two cores consisted essentially of homogenous organic sediment (gyttja) (Fig 3). CON-ISS identified three distinct assemblage zones along the composite stratigraphy. Table 1. Radiocarbon ages obtained from dating of organic sediment from the sediment cores AG1 and AG2. Calibrated ages determined using the IntCal13 calibration curve [33].

Reference
Site and depth (cm) Dated material 14  Zone I (6700-6250 cal yr BP) defines the earliest section of AG1. This relatively short time interval had marked presence of larch (Larix laricina (Du Roi) K.Koch.), balsam fir and, to a lesser extent, black spruce macrofossils. Some minor charcoal peaks were also identified in this zone.
Zone II (6250-2700 cal yr BP) covers a major portion of the core. Charcoal was nearly absent from this zone that was characterized by the predominance of eastern white pine (Pinus strobus (L.)), a more thermophile and less fire-adapted species. In this zone, black spruce, larch, and balsam fir were less common than in Zone I despite some major peaks ca. 3600 and 5800 cal yr BP. Some isolated peaks of jack pine (Pinus banksiana (Lamb.)) were also observed during this period. This zone was also marked by a high number of lepidopteran scale peaks around 4500, 3600, and 2900 cal yr BP. Zone III (2700 cal yr BP to 975 cal yr BP) represents the most recent portion of the lake record. It includes a subzone between 2700 and 2050 cal yr BP representing a transition between zones II and III. In the subzone, every indicator was nearly absent except for lepidopteran scales and some black spruce macrofossils. Between 2050 and 975 cal yr BP, the landscape returned to an early Holocene similar stage characterized by a higher proportion of boreal taxa such as balsam fir, black spruce and larch. Two high magnitude peaks of lepidopteran scales were observed around 1300, 1700 and 2300 cal yr BP. Charcoal, balsam fir, and black spruce macrofossils showed a delay following those lepidopteran peaks. Some remains of jack pine were also identified between 1000 and 2000 cal yr BP, however their abundance was very low relative to the other taxa.
A total of 87 lepidopteran scale peaks and 41 charcoal peaks were identified from both cores (Fig 4). Although not all scales from these peaks could be identified based on the morphotypes described in Navarro et al. [24], 62% of the identified scales corresponded to a SBW morphotype. The ANOVA confirmed that Choristoneura fumiferana was significantly more represented than the other outbreak species (p < 0.0001). Morphotypes of the forest tent caterpillar (Malacosoma disstria (Hübner)) and hemlock looper (Lambdina fiscelaria (Guénée)) were less frequently identified ( Fig 5) and their relative occurence were not significantly different (p = 0.28) from each other. Identification based solely on shape morphotype prevented us  Lambdina fiscelaria occurrences were not significantly different (p = 0.28) from each other. Choristoneura fumiferana was significantly more represented than the other outbreak species (p < 0.0001) and the relative abundance of unidentified scales was significantly higher than the relative abundance of identified scales (p < 0.0001).
https://doi.org/10.1371/journal.pone.0204316.g005 from identifying the majority of those scales (77%) extracted from the peaks as many either were broken or folded.
The frequency of scale peaks was greater than that of charcoal between 5500 and 2500 cal yr BP reaching a frequency of 13 peaks per 1000 years at ca. 3200 cal yr BP. Between 1500 and 1000 cal yr BP there was a reversal of this trend with a fire event rate reaching a maximum frequency of 13 peaks per 1000 years.
The first two axes of the PCA explained 47.6% of the total variation, 32% represented by Axis 1 and 15.6% represented by Axis 2 (Fig 6). The y axis separates Zone 2, mainly influenced by the SBW outbreak frequency and eastern white pine macrofossil abundance, from zones I and III that were more influenced by fire frequency, black spruce and larch macrofossil abundance. Zone I samples were less dispersed than Zone III samples reflecting the higher impact of larch macrofossils. The opposing directions of SBW outbreak frequency and balsam fir macrofossils along with the opposing relationship between fire frequency and eastern white pine macrofossils indicate that these variables varied inversely.

Fire history and climate
The charcoal accumulation rates from Lake Flévy match the known Holocene fire history in eastern Canada [8,[21][22][23]: a more fire-prone early and late Holocene separated by a reduced fire regime-relative to the present-during the mid Holocene. The Holocene fire history recorded in Lake Flévy tracks the δ 18 O-based, southern Ontario postglacial precipitation reconstruction of Edwards et al. [14] as well as different paleofire regime modelisations from eastern Quebec [21,42,43]. Post-and neoglacial periods were cold and dry leading to a higher fire frequency and a dominance of more fire-adapted trees, such as Picea mariana, that are resilient to cooler conditions and drought [44]. Although the mid Holocene was warmer [15], fire activity was very low, likely explained by a high relative humidity during this period [14] and the presence of a more stable air mass leading to less frequent drought events [45].

SBW outbreaks over the Holocene
Based on our dataset, SBW outbreaks are not a rare phenomenon at the Holocene-scale which is in contradiction with our first hypothesis. Multiple peaks of scale abundance were identified downcore. Some have already been previously noted in a mire of the same area (<10km) using faecal macrorests, including outbreaks ca., 1930, 4500, 6650, 7180 and 7560 cal yr BP [9,12,46,47]. However, the novel use of lepidopteran scales also identified multiple, previously unobserved, periods of high lepidopteran abundance such as events at 2500, 2930, 3110, and 3850 cal yr BP.
The scale peak at 4500 cal yr BP has been linked to a possible combined action from SBW and hemlock looper across northern Maine and southern Quebec [9,46,47]. Based on our identification criteria, outbreaks observed in our study site were not associated with hemlock looper. However, given the morphological similarity of the two species' scales, the exceptional quantity of microfossils recovered for this period and the unusual recent hemlock looper activity in the Laurentian Wildlife Reserve (in which Lake Flévy is situated) [48], we cannot completely exclude the possibility of the presence of both defoliators in the area.
The high number of scale peaks identified in the Lake Flévy sediment record allows us to propose the first detailed Holocene-scale SBW outbreak record. The frequency of outbreaks varies over time with periods having a few outbreaks events (10000-6000 cal yr BP; 2500-1000 cal yr BP) and periods marked by very frequent events (5500-2500 cal yr BP).
At a finer temporal scale, some lepidopteran outbreaks (1300, 1800, 4650 cal yr BP) seem to be associated with lagged charcoal, Abies balsamea and Picea mariana macrofossil peaks (Figs  3 and 5). This suggests a direct (and probably combined) influence of these disturbances affecting the input of plant macrofossils into the sediment record.

Fire-insect interactions
Previous studies [49,50] have shown that SBW outbreaks can enhance fire activity over the short term, providing a 5 to 10 yr window of opportunity during which "ladder fuel" builds up due to crown breakage and windthrow of the killed and damaged trees [49,50]. Following this SBW outbreak-related window, humidity rises and the dead wood generated by the outbreak begins to rot [51].
Using a 300 years model of disturbances interactions, James et al. (2011) argue that this ephemeral increase in fire risk due to budworm activity rarely leads to consecutive disturbances at longer timescales [25]. SBW is known to be a cyclic, age-dependent, selective disturbance whereas fire ignition requires a more random event such as a lightning stike. Thus, there is less chance of an immediate succession between these two disturbances. In fact, the two disturbances showed an inverse relationship over the long-term (Fig 6) suggesting an effect of competition for limited resources (mature trees). When fire frequency is low, trees can grow fully and forests can reach a mature stage that enhances SBW activity. On the other hand, when fire frequency increases, there is less time for mature stands to become established, thereby limiting SBW activity.

Conclusion
Our results constitute an important starting point for combining and understanding the longterm interactions between natural disturbances, such as fire and SBW outbreaks, in the boreal forest. This approach provides valuable insight into long-term forest dynamics and may reveal spatial heterogeneity in disturbances (and their interactions). Recognizing these patterns would favor the development of different types of future sustainable management policies depending on the region. For example, this long-term approach could be applied to reconstructing western spruce budworm (Choristoneura occidentalis (Free.)) or mountain pine beetle (Dendroctonus ponderosae (Hopk.)) infestations in western Canada and United States and determining the interactions between outbreaks, fire, and climate in these regions. In a context of global change, it is critical to understand the long-term dynamics of the boreal ecosystem and combined role of multiple disturbances.