Assessing the responses of Sphagnum micro-eukaryotes to climate changes using high throughput sequencing

Study site and field experiment

The study site, Linje peatland (53°11′15″N, 18°18′34″E), is situated in northern Poland within the Complex of Chełmno and Vistula Landscape Parks at 91 m a.s.l. The climate is classified as warm temperate, fully humid with warm summers (Kottek et al., 2006) or cold without dry season with warm summer (Beck et al., 2018). Average annual precipitation in this region is between 500–550 mm and average annual air temperature is between 7.5–8.0 °C. The peatland is dominated by Sphagnum fallax, Eriophorum vaginatum, Oxycoccus palustris, Andromeda polifolia, Aulacomnium palustre, and Drosera rotundifolia (Boiński & Boińska, 2004). The experimental site included also shrubs and young trees such as Betula nana and Pinus sylvestris, respectively.

To simulate warming, we used Open Top Chambers (OTCs), which are commonly used in manipulative experiments in situ (e.g., Delarue et al., 2011; Jassey et al., 2013a; Marion et al., 1997), while water table depth was manipulated by adding or removing peat (Lamentowicz et al., 2016; Reczuga et al., 2018). The advantage of manipulative field experiments is that they are performed under realistic environmental conditions with community composition and abiotic drivers as they appear and work in the field (Underwood, 2009). Field experiments were approved by The Regional Directorate for Environmental Protection in Bydgoszcz (WPN.6205.70.2011.KLD). The experimental design with plots and treatments are given in Fig. 1.

In August 2012, experimental plots were defined according to a full factorial design involving two temperatures and three water table depths. The temperature increase effect was manipulated using OTCs or let under ambient temperature as control (AMB). OTCs are passive warming systems that allow increasing mean temperature by ca. 1.2 °C to 1.8 °C (Marion et al., 1997). The hexagonal OTC chambers were constructed with transparent polycarbonate material. The internal diameter between opposite corners was 250 cm at the base. To allow for air circulation, the OTCs were placed 10 cm above the soil surface. We manipulated water table depth by adding or removing 10 cm layers of peat in 1-m² plots. For each plot four subplots (50 × 50 cm and 30 cm thick) were cut. They were attributed to different measures (e.g., moss sampling for this study, respiration, biomass). Each plot belonged to one of three water table treatments: control (subplots cut and put back in place, abbreviation: CON), wet (10 cm of peat removed beneath each subplot before putting the subplot back in place, abbr.: WET), dry (10 cm of peat added from the wet treatment before putting the subplot back in place, abbr.: DRY). CON treatment was cut and put back in place to control for the disturbance introduced by cutting. Each plot was a combination of one of the three water table treatments and one of the two warming treatments. All six combinations were replicated three times, giving a total of 18 plots (subplots within the experimental plot belonged to the same treatment). Treatments were assigned to the plots randomly. Air temperature was recorded in each plot at 30 cm height using HOBO U23 Pro v2 data loggers (Onset Computer Corporation, USA) and the water table depth (WTD) was measured at each sampling (given in cm below the soil surface). Whenever we are referring to the measured values we use abbreviation WTD, while term “water table treatments” is used when we refer to the levels of manipulation (i.e., DRY, CON, and WET). A subset of the data from the experiment described in this study has already been published elsewhere (Jassey et al., 2018; Lamentowicz et al., 2016; Samson et al., 2018). The experimental design was described by Lamentowicz et al. (2016), where further details can be found. The effect of warming manipulation on temperature has been published in Samson et al. (2018), a study describing influence of warming and water table treatments on carbon dioxide release.

DNA extraction, PCR amplification and next-generation sequencing

Two to three individual Sphagnum stems (top three cm, including capitula) were sampled in three different seasons in each of the plots on May 16th (spring), August 20th (summer) and November 7th (autumn) in year 2013 and fixed in LifeGuard^™ Soil Preservation Solution (MoBio, Carlsbad, CA USA). The plots were placed in monospecific Sphagnum fallax lawns and therefore we did not consider the species level identification of the sampled Sphagnum. In total 54 samples were taken (6 combinations of treatments * 3 replicates * 3 seasons). We sampled only the top part of the mosses as previous studies proved this part to be the most sensitive to the changes in the water table depth and warming manipulations (Basińska et al., 2020; Reczuga et al., 2018). DNA was extracted using the PowerSoil®DNA isolation kit (MoBio, Carlsbad, CA USA) according to the manufacturer’s instruction. The DNA was extracted from two to three sampled Sphagnum mosses, which is approximately 0.3 –0.5 g of material per plot. PCR amplification of the SSU rDNA V9 region was carried out with eukaryotic-specific primers 1380F (CCCTGCCHTTTGTACACAC) and 1510R (CCTTCYGCAGGTTCACCTAC) (Amaral-Zettler et al., 2009). PCR reactions were conducted according to the following conditions: denaturation at 95 °C for 3 min, 25 cycles at 94 °C for 30 s, 57 °C for 30 s and 72 °C for 60 s and final extension at 72 °C for 10 min (Amaral-Zettler et al., 2009). PCR reactions were run with a SensoQuest Labcycler (GmbH, Göttingen, Germany) with 3 l of environmental DNA extract, 4 μL of 5X Colorless GoTaq®Buffer, 0.6 μl of each primer, 0.6 μl of dNTP Mix 10 mM (Promega, Dübendorf, Switzerland) and 0.2 μl of 0.05 U/ μl Go Taq Polymerase (Promega, Dübendorf, Switzerland). Each PCR reaction was performed with a positive and a negative control. Negative controls included PCR reaction mix with DNA extract replaced with water. The purification step was done using the Wizard®SV Gel and PCR Clean-Up System (Promega, Madison, WI USA). Library preparation and paired end sequencing was performed at the iGE3 Genomics Platform of the University of Geneva, Geneva, Switzerland using Illumina Hiseq technology (2 × 150 bp).

Sequence analysis

The computational pipeline used to analyse the reads includes the following steps: trimming of the tagged primer sequences, quality check, removal of rare sequences (i.e., OTU’s occurring less than 3 times in the full dataset), clustering, and taxonomic assignation. For a given read, the quality check was based on moving windows of 50 nucleotides. The probability of incorrect base call was calculated for every nucleotide based on the phred score and arithmetic mean of these probabilities was calculated for every window. To avoid artefactual sequences; for the same purpose, we kept only reads that were present at least three times in the full dataset (De Vargas et al., 2015). Sequence clustering was then performed using Swarm v. 1.2.12 (Mahé et al., 2014). Taxonomic assignation of the resulting OTUs (Operational Taxonomic Units) was done by pairwise alignment of the dominant sequences of each OTU against a selection of dereplicated V9 regions from the PR² database (Guillou et al., 2013) using ggsearch36 v. 36.3.6 (Pearson, 1999). The sequences assigned to Metazoa and Embryophyceae were excluded from further analysis. In order to remove possible false positive sequences only OTU’s occurring at least 10 times in the dataset were retained for further numerical analyses. Based on the taxonomic assignation, each OTU was assigned to a functional group (autotrophs, mixotrophs, parasites, osmotrophs and phagotrophs) according to Table 1. The assignation of OTUs to the functional groups was based on the expert curation (Enrique Lara). OTUs annotated in PR² database with an unknown taxonomic identity (Guillou et al., 2013), were cross-checked individually by aligning them against the NCBI’s nucleotide database using BLAST algorithm with default parameters.

Growing degree day

To link the micro-eukaryotic communities with warming, we used a growing-degree day (GDD) approach, also known as accumulated degree day (ADD). The GDD is widely being used to predict plant and insect phenology and has more recently been used to relate e.g., grasshopper and butterfly community changes to climate change (e.g., Cayton et al., 2015; Nufio et al., 2010). ${\displaystyle \mathrm{GDD}=\frac{T_{\max }+T_{\min }}{2}-T_{\mathrm{base}}}$ where,

T_max are daily maximum air temperature in ^∘C,

T_min are daily minimum air temperature in ^∘C,

T_base is a base temperature, equal to 10 °C.

For each sample, the GDD was calculated for 135 days—the number of days since the beginning of the year until the first sampling.

Numerical analysis

All statistical analyses were performed in R version 3.5.1 (R Core Team, 2012) using RStudio Version 1.1.423 (RStudio Team, 2012). Whenever measured water table depth is used as a variable, we call it WTD, while whenever term “water table treatments” is used we refer to the levels of manipulation (i.e., DRY, CON, WET). The warming treatment effect on GDD and water table treatment on WTD was tested using ANOVA followed by Tukey multiple comparisons of means. We tested for treatment effects (both warming and water table manipulation) on the diversity patterns of micro-eukaryotic communities using the Shannon diversity index, rarefied OTU richness and evenness. OTU richness was estimated using rarefy function in the vegan package (Oksanen et al., 2014) at the sample size of 2410 (minimum number of reads across all samples) to compensate for the variability of the sample sizes. The sequencing effort between the samples was compared by drawing rarefaction curves using rarecurve function and by comparing slope calculated using rareslope function of the vegan R package (Oksanen et al., 2014). The slope was also calculated at the sample size of 2410. To test the significance of differences along measured WTD and temperature gradient (calculated GDD) linear mixed effects models with measured WTD and calculated GDD as fixed factors and season nested in plot as a random effect to account for the differences between plots across time. The assumptions of homoscedasticity and normality were previously tested. The tests were performed using lme function in the nlme package (Pinheiro et al., 2019). To further test the correlation between diversity and WTD and GDD, we calculated Pearson’s correlation using cor.test function.

Subsequently, to test the response of micro-eukaryotes to the GDD, measured water table depth (WTD (cm)) and seasons, at the community level, we performed an RDA (Redundancy Analysis) based on Hellinger-transformed OTU counts. The RDA was done using rda function in the vegan package (Oksanen et al., 2014). Significance tests of each variable and axis in RDA models were tested using permutation tests (anova.cca function in the vegan package). To identify micro-eukaryotic indicators of climate change, we used Dufrene-Legendre indicator species analysis (Dufrêne & Legendre, 1997) and multipatt function in the indicspecies package (De Cáceres & Legendre, 2009). To manually verify taxonomic assignation against more wholesome database, the selected indicators’ (with the highest indicator value) sequences were cross-checked by aligning them against the NCBI’s nucleotide database using BLAST algorithm with default parameters.

Warming and water table manipulation effects on microclimatic conditions

Mean maximum daily temperature was about 1.1–1.2 °C higher in the OTC than in the ambient plots (Samson et al., 2018). The warming effect varied seasonally and over daytime, reaching a maximum effect of Δ 4.1 °C on 5th of May 2013 and causing a slight cooling effect during winter and during the night (Lamentowicz et al., 2016). Mean daily air temperature was calculated on the basis of 10-minute averages. For each sampling, the temperature of the month prior to the sampling was calculated as the average of the mean daily air temperature. This mean was higher in the OTC than in ambient plots in spring and summer periods (by 0.54° and 0.43°, respectively). In each season WTD was higher (i.e., lower water table) in DRY as compared to the CON and WET. Mean WTD in DRY was 15 cm, 46 cm and 25 cm in spring, summer and autumn sampling, respectively, while in CON it was nine cm, 40 cm and 19 cm and in WET it was eight cm, 40 cm, 19 cm (Table S1, Fig. S1). Significant differences were observed between CON and DRY as well as between WET and DRY (Tukey, p < 0.05), while no significant difference was observed between WET and CON (Tukey, p = 0.6). GDD differed significantly (ANOVA, p <0.05) between AMB and OTC across seasons, with mean GDD of 76, 732 and 549 in AMB spring, summer and autumn periods, respectively, while OTC mean GDD was 101, 847 and 615 in spring, summer and autumn periods, respectively (Table S2). Significant differences between GDD in AMB and OTC were observed at each season (Tukey, p < 0.05). The environmental variables in experimental plots are provided in Data S1.

Micro-eukaryotic community structure

The Illumina sequencing generated 37,804,202 raw eukaryotic SSU V9 reads. Of these, the number dropped to 29,791,631 after quality check. These sequences were then clustered in 2,043 OTUs, reduced to 1,397 OTUs after removing sequences assigned to Metazoa or Embryophyceae. The abundance of each of the 2,043 OTUs in experimental plots is shown in Data S2, while the dominant sequence of each OTU is shown in Data S3. For ecological analysis, all assignations with similarity below 80% were excluded, leaving 1,226 OTUs. The majority of OTUs were assigned to the following five eukaryotic supergroups in decreasing order of taxonomic richness (i.e., 614, 214, 152, 99 and 83 OTUs): Opisthokonta, Alveolata, Archaeplastida, Rhizaria, respectively. A total of 369 OTUs (26.4% of the total number of OTUs) were assigned to taxa with a similarity above 97%. Among these, 250 were assigned to fungi and 119 to other taxa. Fungi OTUs reaching 97% of identity with the database constituted 40% of all fungi while only 16% of the other taxa reached the identity threshold illustrating the knowledge gap between fungi and protists.

Warming and water table manipulation effects on micro-eukaryotic diversity

The rarefaction curves (Fig. S2) clearly showed that the sequencing effort was not saturated at the sample level. However, when considering the whole dataset, the rarefaction curve indicated that the majority of micro-eukaryotic diversity was recovered. To compensate for the saturation variability among treatments, the OTU richness was calculated on a rarefied community in order to compare the effect of the different treatments. We did not observe significant differences in the slope of the rarefaction curves among water table or warming treatments (ANOVA, p = 0.6 and p = 0.9, for warming and water table treatments, respectively; Fig. S3). Linear mixed effect models revealed that Shannon diversity, evenness and rarefied OTU richness differed significantly along WTD and GDD (ANOVA, p < 0.01, Table S3). All diversity indices declined significantly with increasing WTD (Fig. 2, cor = −0.3 for all diversity measures, p < 0.05). No significant correlation was observed between GDD and biodiversity metrics. Shannon diversity and OTU richness were higher in autumn as compared to summer and spring. No significant difference was observed between warming treatments. The Shannon diversity, evenness, OTU richness and rarefied OTU richness in each analysed sample are given in the Table S4.

Micro-eukaryotic community response to the experimental manipulation

The redundancy analysis (RDA) revealed significant correlations between micro-eukaryotic communities and WTD (anova.cca permutation test, p = 0.001), GDD (p = 0.001) and season (p = 0.001). The full RDA model explained 22.5%, most of this (16.6%) being explained by seasons. To disentangle the seasonal effects from the influence of other environmental factors, we performed RDA with season as a conditional variable (Fig. 3). Here again, WTD (p = 0.014) and GDD (p = 0.017) were significant. The partial RDA model explained 5.9% (p = 0.003) of the variance. RDA model with water level and temperature treatments as variables and with season as a conditional variable revealed significant differences between water level treatments (p = 0.001), while the difference between temperature treatments was insignificant (p = 0.147; Fig. S4).

Among the 1,226 OTUs, 216 were considered autotrophs, 52 mixotrophs, 42 parasites, 579 osmotrophs and 319 phagotrophs (Table 1). Eighteen OTUs that could not be placed phylogenetically within any phylum, and therefore could not be functionally assigned, were excluded from the analysis. In the DRY treatment in spring, we observed a higher relative abundance of autotrophs and phagotrophs and a lower relative abundance of osmotrophs and parasites as compared to the WET and CON treatments. No such difference between treatments was observed in summer and autumn (Fig. 4).

Micro-eukaryotic indicators of global change

The indicator species analysis with each treatment as a separate group, allowed identifying 109 indicator OTUs (Table S5). A total of 96 OTUs were DRY indicators, 10 were WET indicators and three were identified as CON indicators (p < 0.05). Among the 96 DRY indicators, 28 were assigned to Fungi, 23 to Chlorophyta, 17 to Ciliophora, 14 to Cercozoa, six to Dinophyta, two to Cryptophyta and Stramenopiles, one to Apicomplexa, Discoba, Mesomycetozoa and Streptophyta. The DRY indicator with the highest indicator value was OTU_X77 (100% similarity with Cortinarius sp.) which had an indicator value of 0.963. The relative abundance of this OTU was highest in autumn (Fig. 5). Among the 10 DRY indicators with the highest indicator value, five were assigned to Archaeplastida (including four assigned to Trebouxiophyceae), four to Fungi and one to Alveolata (Colpodidae; Table 2). The increase of the relative abundance of autotrophs in spring corresponded to indicator algae adapted to dry soil environment and capable of symbiosis with lichens, such as Elliptochloris, (Trebouxiophyceae, Metz et al., 2019). The relative abundance of OTUs assigned to Elliptochloris, predominantly terrestrial taxa (Gustavs et al., 2017; Rindi, Hodkinson & Jones, 2011), increased from 1% in CON plots in May to 6% in DRY plots in May (Fig. S5), showing a shift towards a more terrestrial community.

A total of 27 indicators were identified for ambient vs OTC, nine of which were indicators of ambient temperature and 18 were indicators of warming (p < 0.05; Table S6).