Niche Conservatism Drives the Elevational Diversity Gradient in Major Groups of Free-Living Soil Unicellular Eukaryotes

Ancestral adaptations to tropical-like climates drive most multicellular biogeography and macroecology. Observational studies suggest that this niche conservatism could also be shaping unicellular biogeography and macroecology, although evidence is limited to Acidobacteria and testate amoebae. We tracked the phylogenetic signal of this niche conservatism in far related and functionally contrasted groups of common soil protists (Bacillariophyta, Cercomonadida, Ciliophora, Euglyphida and Kinetoplastida) along a humid but increasingly cold elevational gradient in Switzerland. Protist diversity decreased, and the size of the geographic ranges of taxa increased with elevation and associated decreasing temperature (climate), which is consistent with a macroecological pattern known as the Rapoport effect. Bacillariophyta exhibited phylogenetically overdispersed communities assembled by competitive exclusion of closely related taxa with shared (conserved) niches. By contrast, Cercomonadida, Ciliophora, Euglyphida and Kinetoplastida exhibited phylogenetically clustered communities assembled by habitat filtering, revealing the coexistence of closely related taxa with shared (conserved) adaptations to cope with the humid but temperate to cold climate of the study site. Phylobetadiversity revealed that soil protists exhibit a strong phylogenetic turnover among elevational sites, suggesting that most taxa have evolutionary constraints that prevent them from colonizing the colder and higher sites of the elevation gradient. Our results suggest that evolutionary constraints determine how soil protists colonize climates departing from warm and humid conditions. We posit that these evolutionary constraints are linked to an ancestral adaptation to tropical-like climates, which limits their survival in exceedingly cold sites. This niche conservatism possibly drives their biogeography and macroecology along latitudinal and altitudinal climatic gradients.


Introduction
Niche conservatism is a hypothesis proposed to explain why biological entities exhibit diversity (biogeographic) patterns along climatic gradients [1]. Niche conservatism states that taxa remain in the same climate in which their ancestors evolved because they tend to retain their ancestral niche traits over the evolutionary course [1,2]. This would explain why most taxa fail to disperse into new climates, providing an underlying mechanism for the biogeographic patterns observed along latitudinal or altitudinal climatic gradients [2,3].
Since most modern taxa originated in tropical-like climates [4], niche conservatism predicts that biodiversity will peak in warm, humid (optimal) climates because they should have evolutionary constraints that prevent them from thriving at exceedingly hot or cold (suboptimal) climates [1]. Accordingly, biodiversity will be better predicted by air temperature at sites with humid and cold climates and by water availability at sites with hot and dry climates (Prediction 1).
Niche conservatism also predicts that few taxa have evolved new traits (evolutionary novelties) to cope with suboptimal climates [1,5]. If so, most taxa should have narrow geographic ranges restricted to optimal climates, while few taxa should have broad geographic ranges extending beyond optimal climates [6]. This idea is consistent with the Rapoport effect [6,7], a macroecological pattern that describes an increase in taxa ranges from optimal to suboptimal climates (Prediction 2).
Due to niche conservatism, closely related taxa share ancestral niche traits (Prediction 3). Therefore, these taxa cannot coexist in places with optimal climates because they compete for similar resources [5,8]. Optimal climates then promote community assembly via competitive exclusion, resulting in phylogenetic overdispersion or in the coexistence of distantly related taxa with non-overlapping niches [8,9]. However, closely related taxa that evolved and share traits to survive in cold or hot conditions can coexist in suboptimal climates [8]. This occurs because suboptimal climates filter out taxa that are not adapted to survive in cold or hot conditions [10]. Therefore, closely related taxa adapted to suboptimal climates compete less with each other under cold or hot conditions because they have access to more resources than in optimal climates [8,10]. Suboptimal climates then promote community assembly via habitat filtering, resulting in phylogenetic clustering or in the coexistence of closely related taxa sharing evolutionary novelties to cope with exceedingly hot or cold climates [8,9].
Finally, since evolutionary constraints prevent taxa from colonizing suboptimal climates [1], habitat filtering will limit the dispersal of less adapted taxa as they approach the hotter or colder end of the climatic gradient [11,12]. Phylogenetic beta diversity (phylobetadiversity) or phylogenetic variation among communities [13] should therefore reveal changes in phylogenetic composition along any climatic gradient (Prediction 4).
While these predictions are consistent with niche conservatism, they are also consistent with the idea that there are abiotic constraints on the conditions that an organism can readily tolerate (i.e. they do not require that these be constrained by the past evolutionary history of the group). Therefore, local physicochemical conditions represent a confounding variable that needs to be controlled for in order to examine a cause-and-effect relationship between niche conservatism and the occurrence of diversity patterns in nature [14,15]. Typically, this is achieved by using a restrictive sampling strategy, which consists of assessing predictions along a climatic gradient of comparable habitats to reduce the confounding effects of local abiotic factors [9,14,15].
Niche conservatism predictions have been extensively tested in multicellular organisms [1][2][3][4], showing the usefulness of niche conservatism in explaining biogeographic patterns. By contrast, they have been seldomly tested in unicellular organisms, including bacteria [5,9,10] and protists such as testate amoebae [12]. Thus, many biogeographic generalizations proposed for plants and animals have not yet been sufficiently tested for microorganisms, particularly for soil protists [16].
Soil protists are good model organisms for testing whether niche conservatism applies to unicellular organisms. Soil protists appear to retain their traits over their evolutionary course, including a need for warm, humid climates [12,17]. Evidence suggests that at least some soil protist taxa evolved under the tropical-like climate that characterized much of the Mesozoic [18,19]. These past climatic conditions seem to have imposed evolutionary constraints on their ability to adapt to exceedingly hot or cold climates, influencing their current diversity patterns, performance and fitness [12]. Indeed, soil protists diversity peaks in warm, humid climates [15,20,21] and declines towards exceedingly hot [22] and cold [23] climates. They also show high survival, growth and reproductive rates at warm temperatures (ca. 18 °C and 24 °C) than at exceedingly cold or hot temperatures [24][25][26]. This background suggests that most soil protists lack ecophysiological adaptations to survive and reproduce in exceedingly cold or hot climates. Probably, the retention of their ancestral adaptation to tropical-like climates imposes strong constraints on their ability to colonize sites with suboptimal climates and drives their current biogeographic patterns.
Herein, we tested niche conservatism predictions on free-living representatives of five major, phylogenetically and functionally contrasted common soil protist groups: Bacillariophyta (diatoms), Cercomonadida (cercomonadids), Ciliophora (ciliates), Euglyphida (euglyphids or filose testate amoebae) and Kinetoplastida (kinetoplastids) and a group combining all these taxa (herein refer to as the "soil protist group"). We selected these soil protist taxa because (1) they are frequent, abundant and diverse in soils; (2) they have different life history strategies and functional roles; and (3) they are well distributed across the eukaryotic tree [27]. Our rationale behind this choice was to be able to generalize conclusions to as many soil eukaryotes as possible.
Niche conservatism predictions were tested along an elevational gradient in Switzerland, which exhibits moderate to high humidity over the year but decreasing temperature with altitude [28]. This elevational gradient is ideal for assessing niche conservatism in soil protists because the forests are dominated by Beech (Fagus sylvatica) throughout the gradient [29]. Beech trees require moderate to high precipitation and through shading and litterfall establish similar microclimate and soil properties on sites where they are the dominant plant species [30,31]. Therefore, beech-dominated forests represent comparable habitats for unicellular organisms [31]. This feature reduces the confounding effects of local abiotic factors [14,15] and contributes to unravelling the role that niche conservatism plays in shaping soil protist diversity patterns.

Study Site and Sampling Strategy
The study was conducted in Western Switzerland, which has a temperate maritime climate [28] (Fig. 1). Within this region, forests are mostly dominated by beech (F. sylvatica) trees, and, although mostly exploited for timber, their structure and vegetation composition remains similar to pristine forests [29]. The studied sites were selected within the permanent plots of the Swiss Biodiversity Monitoring (BDM) program (https:// bit. ly/ 3puxE xC). We chose 10 plots from 458 to 1308 m a.s.l. to study biodiversity along nearly the entire elevation range of beech-dominated forests in Switzerland. This elevational gradient exhibits high humidity over the year and decreasing temperature with elevation (mean annual precipitation and temperature: 1200 mm and 8.3 °C, respectively) [28,29].
Beech-dominated forests have comparable microclimate/ soil properties and constitute similar habitats for microorganisms [30,31]. We therefore collected samples at the periphery of the BDM plots ( Fig. 1) to reduce the confounding effects of abiotic factors [9,14,15]. At each site we collected three soil cores (5-cm diameter × 5-cm depth) and then mixed them to obtain a composite sample for each sampling site. We kept the soil samples cool at all times and performed DNA extraction within 24 h of sampling.

Characterization of Soil Protist Communities
We followed standard protocols published elsewhere [32] for DNA extraction, PCR and sequencing (Illumina, targeting the SSU rRNA gene V9 region of eukaryotes), as well as for taxonomic assignment of the reads obtained. These protocols are summarized in the supplementary material ESM1.
We randomly subsampled 50,000 reads from each composite sample to account for unequal sample sizes and defined a community as all operational taxonomic units (OTUs) originating from a single sampling site. From these subsamples, we selected those assigned to Bacillariophyta, Cercomonadida, Ciliophora, Euglyphida and Kinetoplastida. Another group (the soil protist group) was also created combining all the reads assigned to those taxa. We estimated taxonomic richness in each community based on the total number of OTUs, a metric that in turn was correlated to phylogenetic diversity (see ESM1: Fig. S1). We assessed the quality of our sampling effort by building species accumulation curves for each taxon in the package vegan [33] and using R 3.1.2 [34].

Abiotic Factors and Confounding Effects
Abiotic factors (soil temperature, moisture, humus content, continentality, light, nutrients, pH) were estimated at each sampling site using Landolt indicator values [35]. This bioindication method estimates soil abiotic factors by averaging the individual indicator values (i.e. the individual ecological tolerances) of plant species present on a site. The average indicator values range between 1 and 5 and provide robust information on the long-term environmental conditions characterizing a site [36]. Abiotic factors were estimated using Landolt indicator values because their performance in predicting biodiversity patterns outperforms that of data interpolated from climatic databases popularly used in ecological studies [37]. In part, this is because climate databases still exhibit data gaps and low resolution, particularly in forested regions of Europe such as our study site [38].
After reducing the set of abiotic factors to a representative subset to avoid multicollinearity (we retained soil temperature, light, humidity, pH and humus content), we constructed linear models with protist diversity as response variable and the subset of abiotic factors as descriptors and tested the overall model significance as well as the significance of each descriptor using ANOVA in the package vegan [33]. Since our sampling strategy aimed to decrease the confounding effects of abiotic factors to disentangle the role of niche conservatism in shaping protist diversity patterns, we expected to record a lack of correlation between abiotic factors and soil protist diversity.

Prediction 1
The elevational gradient investigated exhibits high humidity but decreasing air temperature with altitude [28,29]. Therefore, according to niche conservatism, temperature should limit protist survival along this humid but cold altitudinal gradient. To test this prediction, we standardized OTU richness to avoid the unwanted effects of sampling artefacts. This process was conducted using a combination of range interpolation and species richness estimates. Range interpolation consisted in dividing the entire elevational gradient into 100-m elevational bands and assuming that taxa were present at all elevations between the lowest and highest observed elevations [39]. Richness estimates involved the calculation of the expected OTU richness for each 100-m elevational band using the richness estimator available in the package iNEXT [40]. Then, we estimated the mean annual temperature for each 100-m elevational band based on the mean annual air temperature reported for the lowest sampling site (10 °C) [41] and assuming a decrease (lapse rate) of 0.6 °C for each 100-m increase in elevation [39]. While the estimated and the empirical rates at which air temperature cools with altitude may differ, the use of a lapse rate may not pose an issue in our case. Our estimate of air temperature decreases monotonically with altitude, mirroring the empirical trend observed at the study site in a previous study [28]. We then constructed linear and quadratic models using richness as response variable and both the elevation and the mean annual air temperature as descriptor variables. The best-fitting model was selected based on the Akaike's information criterion for each model.

Prediction 2
We used the midpoint method [6] to assess the elevational Rapoport effect. This method analyses the correlation between elevation and the mean elevational range size of all taxa present at each sampling site. The correlation between both variables was tested using Pearson's correlations in R 3.1.2 [34]. A positive correlation between both variables would support the occurrence of an elevational Rapoport effect, i.e. an increase in the elevational range sizes of soil protists from lower and warmer sites to higher and colder sites. This outcome would suggest that soil protists are evolutionarily constrained to cope with cold climates.

Prediction 3
Analysis of phylogenetic structure informs whether taxa coexisting within a community are distantly related (or phylogenetically overdispersed) or closely related (or phylogenetically clustered) [8]. We therefore investigated phylogenetic structure within communities to assess the relatedness among taxa coexisting within the same community. For this task, we used two indices, namely, the mean pairwise distance index (-NRI) and the mean nearest taxon distance index (-NTI) [42]. -NRI estimates the mean phylogenetic distance among all pairs of taxa within a community, whereas -NTI estimates the mean phylogenetic distance between each taxon and its closest relative within a community. Positive values of these indices indicate phylogenetic overdispersion, and negative values indicate phylogenetic clustering within communities. If climate is optimal for soil protists along the elevational gradient, communities will be phylogenetically overdispersed [8,9]. Conversely, if climate is suboptimal for soil protists, then communities will be phylogenetically clustered [8,9]. Finally, if soil protists do not conserve their traits over evolution, their within-community phylogenetic structure will be predicted by chance [1,42]. Both -NRI and -NRI indices were estimated using the package picante [42], and their outcomes were compared against values predicted by a null model computed using a trial swap algorithm (999,000 randomizations).

Prediction 4
We investigated phylobetadiversity to assess the elevational phylogenetic variation among communities in the study site. Phylobetadiversity was estimated using the PhyloSor index, which estimates the phylogenetic composition shared by two or more sites [9]. We also decomposed phylobetadiversity in its two additive components, namely, phylogenetic turnover and phylogenetic nestedness [13]. Measuring both components allows defining the proportion of phylobetadiversity that is caused by phylogenetic replacement among sites (or true turnover) and the proportion that is caused by loss of phylogenetic diversity among sites, respectively. If soil protists exhibit a phylogenetically conserved response to climate, then the elevational decrease in temperature will progressively limit soil protist colonization towards colder and higher elevations. If so, phylogenetic composition will vary among communities, and phylobetadiversity will deviate from that expected by chance. In turn, if soil protists do not conserve their ancestral niche traits, then their phylogenetic composition will vary randomly along the elevational gradient since all taxa will be equally likely to colonize the entire gradient of comparable (beechdominated) habitats. The PhyloSor index and its two additive components were estimated in betapart [43], and their outcomes were compared against values predicted by a null model computed using an independent swap algorithm (50,000 randomizations).

Characterization of Soil Protist Communities
We recorded a total of 1413 OTUs and 47,394 reads along the elevational gradient. The Cercomonadida (373 OTUs: 12,437 reads) and Ciliophora (582 OTUs: 15,599 reads) were the most diverse and abundant heterotrophic taxa at each forest site, followed by Kinetoplastida (223 OTUs: 10,377 reads) and Euglyphida (157 OTUs: 8058 reads). The phototrophic taxon Bacillariophyta was the less diverse and abundant (78 OTUs: 923 reads). Species accumulation curves confirmed that sampling effort and sequencing covering were enough to record a significant proportion of the total OTU richness found at each sampling site (supplementary material ESM2: Fig. S1).

Abiotic Factors and Confounding Effects
Soil abiotic factors as inferred from Landolt indicator values varied within narrow ranges among beech-dominated forests (ESM2 : Table S1). An ANOVA test revealed that abiotic factors were not significantly correlated to soil protist diversity ( Table 1), suggesting that our restrictive sampling strategy succeeded in reducing the confounding effects of local abiotic factors.

Prediction 1
Protist diversity exhibited a negative correlation with elevation and a positive correlation with climate (mean annual air temperature) ( Table 2). Therefore, protist diversity peaked at lower and warmer elevations. This was consistent with Table 1 Correlations between soil protist diversity and local abiotic factors estimated at each sampling (beech forest) site AIC Akaike information criterion, F moisture, H humus content, L light, R soil reaction (pH), T temperature After reducing the set of local abiotic factors to a representative subset to avoid multicollinearity (we retained soil temperature, light, humidity, pH and humus content), we constructed linear models with richness as response and variables as descriptors and tested the overall model significance as well as the significance of each descriptor using ANOVA in the R package vegan [33]. P values are shown for the whole model and for each abiotic factor. The abiotic factors were coded (i.e. F, H, L, R,T) using the codes originally proposed by Landolt [35]. The codes are explained in the lower margin of this the idea that most soil protists are subjected to evolutionary constraints which make it challenging for them to cope with colder climates.

Prediction 2
We recorded a positive and significant correlation between elevation and the mean elevational range size of soil protists (Fig. 2). Thus, the elevational range of soil protists increased in size from lower and warmer to higher and colder elevations. This increasing pattern in the elevational range size of soil protists is consistent with an elevational Rapoport effect and suggests that most soil protists have ecophysiological constraints to adapt to cold climates.

Prediction 3
-NRI and -NTI indices varied in their absolute values (ESM2 : Table S2) but led to the same conclusions in terms of phylogenetic structure within soil protist communities (Fig. 3). The analysis of the phylogenetic structure within the "soil protist group" combining all investigated taxa revealed that overall, soil protists exhibit both phylogenetically overdispersed communities composed of distantly related taxa

Prediction 4
Phylobetadiversity as measured by the PhyloSor index revealed that all taxa exhibited a significant dissimilarity in phylogenetic composition between any pair of soil protist communities taken at random from our study site (PhyloSor ranged from 0.57 to 0.66, P < 0.05; Table 3). Phylobetadiversity was mainly the result of phylogenetic turnover (PhyloSor Turn ranged from 0.5 to 0.60, P < 0.05; Table 3) and, to a lesser extent, phylogenetic nestedness (PhyloSor PD ranged from 0.04 to 0.09, P < 0.05; Table 3).

Discussion
Consistent with Prediction 1, soil protist diversity was significantly correlated with climate (mean annual air temperature). This so-called species-energy effect [44] was previously shown to be as an important driver of soil protist biogeography in latitudinal and elevational gradients exhibiting humid climates and a significant thermocline [12,15,21,23]. Other studies have also found that soil protist diversity patterns result from the synergy between climate and local abiotic factors [45,46]. Indeed, this synergy is important in shaping the biogeography of any taxon, although it is also known that the relative importance of climate and local abiotic factors varies significantly with the spatial scale of the study [47]. One way to disentangle the role of climate as a driver of biodiversity is to use a restrictive sampling strategy to reduce the confounding effects of local factors [9,14,15]. In our case, we restricted our sampling to beech-dominated  [30,31]. This restrictive sampling strategy proved to be effective in our case, as soil protist diversity did not exhibit any correlation with local abiotic factors. Instead, the elevational variation of soil protist diversity was significantly correlated to climate (mean annual air temperature). The evaluation of Prediction 2 supported the idea that climate is a good predictor of the elevational diversity gradient of soil protists. All taxa exhibited an elevational Rapoport effect or a progressive increase in their distribution ranges from lower and warmer sites to higher and colder sites. The Rapoport effect has been related to the existence of evolutionary constraints that prevent multicellular organisms [6,7,11] and soil protists [12] from adapting to and colonizing areas with severe, often exceedingly cold/hot climates. So, the observed elevational Rapoport effect suggests that overall, soil protists lack the adaptations needed to cope with the cold climate of higher elevation.
Following Prediction 3, the assessment of within-community phylogenetic structure in the soil protist group revealed both phylogenetically overdispersed communities assembled Fig. 2 Rapoport effect or positive correlation between elevation and mean elevational range size (MER) of all OTUs we found co-occurring at each sampling site (circles). The positive correlation observed between both variables shows that all eukaryotic microbial groups exhibited a significant increase in the size of their elevational ranges with elevation by competitive exclusion and phylogenetically clustered communities assembled by habitat filtering. Since phylogenetic overdispersal and phylogenetic clustering are respectively expected to occur in optimal and suboptimal climates [1,2], we can conclude that the humid, cold climate prevailing at our study sites represents a nearly optimal climate for some taxa and a harsh climate for others. Soil protist taxa exhibiting phylogenetically overdispersed communities probably have a broader thermal tolerance than those exhibiting phylogenetically clustered communities. This outcome is consistent with niche conservatism [1][2][3][4][5] and suggests that soil protists conserve their niche traits among closely related taxa and over evolutionary time.
Assessment of within-community phylogenetic structure in monophyletic taxa revealed that the autotrophic taxon Bacillariophyta (diatoms) exhibited phylogenetically overdispersed communities assembled by competitive exclusion, while heterotrophic taxa Cercomonadida, Ciliophora, Euglyphida and Kinetoplastida exhibited phylogenetically clustered communities assembled by habitat filtering along Fig. 3 Within-community phylogenetic structure as measured by -NRI (circles) and -NTI (squares) indices at each sampling site. Positive values (above zero) indicate phylogenetic overdispersion or communities of distantly related taxa assembled by competitive exclusion. Negative values (below zero) indicate phylogenetic clustering or communities of closely related taxa assembled by habitat filtering. Within-community phylogenetic structure estimations were significant in all cases (P < 0.05) the humid and cold elevational gradient of beech forests. In these forests, diatoms assembled their communities by competitive exclusion because they probably compete for the scarce sunlight available. In fact, sunlit spots are often sparse and scattered in beech forests because the trees develop a dense canopy that limits the amount of sunlight reaching the ground [29]. Therefore, closely related diatoms with similar niche (light intensity/wavelength) requirements must actively compete for this limiting resource. On the other hand, heterotrophic protists possibly assemble their communities by habitat filtering because the prevailing cold climate restricts most of their representatives. Evidence suggests that heterotrophic protists have lower fitness and performance (reproductive, growth and survival rates) than autotrophic protists when exposed to cold temperatures [24][25][26]. Some do not survive or cannot reproduce below 10 °C [24]. Thus, our study site must harbour only a subset of the total pool of heterotrophic taxa found in the region. These closely related taxa must share evolutionary novelties that allow them to cope with the elevational decrease in temperature that characterizes the study site.
Evaluation of phylobetadiversity (Prediction 4) revealed that the phylogenetic composition varies significantly between forests (as shown by PhyloSor). In particular, most OTUs occurred at only one site (as shown by PhyloSor Turn ), while a comparatively minor part occurred at two or more sites (as shown by PhyloSor PD ) [13]. Thus, there is a strong phylogenetic turnover between sites [43], suggesting that each taxon (i.e. Bacillariophyta, Cercomonadida, Ciliophora, Euglyphida and Kinetoplastida) is mostly represented by individuals adapted to survive in specific temperature ranges. Therefore, most soil protists exhibit evolutionary constraints that do not allow them to colonize and thrive in sites that exhibit temperatures that exceed their tolerance range. This lends further support to the existence of evolutionary constraints that prevent taxa from adapting to and surviving in new climatic conditions [11,12]. Thus, niche conservatism represents a potential underlying mechanism for soil protist phylobetadiversity. These results suggest that niche conservatism not only acts at the community level (as revealed by analyses of within-community phylogenetic structure), but also at the metacommunity level.

Conclusion
Previous research and the present study support the idea that soil protist diversity increases towards warm (and humid) climates. However, the novelty of our study lies in the fact that we provide evidence suggesting that this species-energy effect may have an evolutionary origin rooted in the conservatism of ancestral thermal regimes rather than on present-day climatic conditions. The phylogenetic signal of this thermal niche conservatism was herein explicitly detected by assessing the within-community phylogenetic structure (phylogenetic overdispersion/clustering) and phylogenetic beta diversity (phylobetadiversity) along an elevational and climatic gradient of comparable habitats. The outcomes provided by these metrics suggest that soil protists exhibit evolutionary constraints that do not allow them to adapt to climates departing from warm and humid conditions. This thermal constraint might also be driving their current biogeographical and macroecological patterns on Earth and might be the reason why temperature often arises as an important predictor of soil protist diversity over latitudinal and elevational gradients. Considering that the phylogenetic niche conservatism has contributed to successfully explain the occurrence of diversity patterns in plants and animals, we could also state that our study contributes with additional evidence to demonstrate that eukaryotic multicellular and unicellular diversity patterns might be produced and maintained by similar processes. Our findings contribute thus to generalizing broad evolutionary mechanisms to the whole domain Eukarya and, arguably, to all life on Earth. Variation in phylogenetic composition among sites (PhyloSor) and its underlying phenomena, namely phylogenetic turnover (PhyloSor Turn ) and phylogenetic nestedness (PhyloSor PD ), were estimated for each soil protist taxa using betapart [43]. All taxa exhibited higher than expected phylobetadiversity, primarily represented by phylogenetic turnover among sites