Climate warming promotes species diversity, but with greater taxonomic redundancy, in complex environments.

Climate warming is predicted to alter species interactions, which could potentially lead to extinction events. However, there is an ongoing debate whether the effects of warming on biodiversity may be moderated by biodiversity itself. We tested warming effects on soil nematodes, one of the most diverse and abundant metazoans in terrestrial ecosystems, along a gradient of environmental complexity created by a gradient of plant species richness. Warming increased nematode species diversity in complex (16-species mixtures) plant communities (by ~36%) but decreased it in simple (monocultures) plant communities (by ~39%) compared to ambient temperature. Further, warming led to higher levels of taxonomic relatedness in nematode communities across all levels of plant species richness. Our results highlight both the need for maintaining species-rich plant communities to help offset detrimental warming effects and the inability of species-rich plant communities to maintain nematode taxonomic distinctness when warming occur.

INTRODUCTION

Increases in Earth’s surface and air temperatures are predicted to alter the interactions among species (). The alteration of species interactions (for example, predation or competition) is due mainly to differences in temperature-dependent biological rates, such as growth and reproduction, among interacting species (). Species that are able to meet their enhanced biological requirements, such as greater resource consumption due to elevated metabolic demands () at elevated temperature, are likely to become competitively superior to those that are unable to meet enhanced resource requirements (, ). As a consequence, warming may trigger extinction events at local scales, particularly in resource-limited environments (, , ).

Diversity of primary producers or resource species (for example, plants or algae) has been consistently shown to increase primary production in ecosystems (, ). Greater productivity associated with a higher number of resource species often benefits their consumer species, both in terms of species abundances and diversity (–). Haddad et al. () revealed multiple mechanisms causing positive resource diversity–consumer diversity relationships, such as different consumers specializing on a variety of resources or greater resource biomass supporting a greater number of consumer individuals and thereby increasing consumer diversity. A recent theoretical study also revealed that higher availability of resources could dampen warming-induced competitive hierarchies among consumer species and rescue species from local extinction (). However, empirical evidence of these effects is scarce.

Climate warming was recently shown to increase primary production the most in the most biodiverse plant communities (plant community with 16 species), but negligibly so in plant monocultures in a grassland experiment (). These potential interactive effects of plant species richness and warming on primary production may have cascading effects on community structure at higher trophic levels due to alterations in bottom-up forces. Alongside, warming may also directly weaken top-down forces by reducing the population size of predators (). The loss of predators could increase the population sizes of dominant prey species within the community, potentially leading to competitive exclusion of rare prey species (, ). Because communities are often structured by the combined effects of top-down and bottom-up forces (, –), and warming simultaneously affects both these forces (), consumer communities may respond to warming in unexpected ways (, –). For instance, a recent experimental study reported positive effects of warming on species diversity of phytoplankton species and suggested greater top-down regulation of dominant phytoplankton prey species in warm and productive environments as one of the underlying mechanisms ().

Here, we investigate the responses of free-living soil nematode species to experimental warming along an experimental gradient of plant species richness as a representation of simple-to-complex soil environments (). Free-living nematodes are highly abundant organisms in the soil and one of the most diverse metazoan taxa in terrestrial ecosystems (). Soil nematodes heavily depend on plant productivity (), and high plant diversity has been shown to increase the quantity and diversity of resources for soil nematodes via increasing plant and microbial biomass production and rhizodeposition (, –). An important feature of nematode communities is that they are composed of major trophic groups, such as herbivores, detritivores, and predators (). This makes nematodes a unique model taxon to study how warming effects may vary among different trophic groups within a community (, ).

We hypothesized that negative effects of warming on nematode diversity would be greatest in plant monocultures (simpler soil environment) and progressively smaller in more diverse plant communities (more complex soil environments). The basis of this hypothesis is that the more complex soil environments of diverse plant communities would allow more competing nematode species to coexist via greater resource diversity and availability () and thus relax warming-induced competitive hierarchy. In contrast, warming of a simpler soil environment would promote competitive displacement within nematode communities due to lower resource availability and a more homogenous environment (, ). In addition, we tested whether warming- and plant species richness–induced changes in nematode diversity are associated to changes in potential bottom-up (plant biomass) and/or top-down regulation (predator density) of nematode communities.

RESULTS

Plant species richness and warming effects on soil nematodes

Our analyses show a significant interaction effect of plant species richness and experimental warming on nematode taxa richness [coefficient (β) = 0.02, t = 2.47, P = 0.01], that is, nematode taxa richness increased (~38%) in warmer mixed plant communities, but declined (~28%) in warmer plant monocultures (Fig. 1A). Effective number of nematode species showed the same pattern (~36% increase in warmed species-rich plant communities and ~39% decrease in warmed plant monocultures) as nematode richness (β = 0.03, t = 3.24, P < 0.01; Fig. 1B). In contrast, Pielou’s evenness decreased with experimental warming independent of plant species richness (β = −0.20, t = −2.17, P = 0.03; Fig. 1C). Rarefied nematode richness also decreased with experimental warming independent of plant species richness (β = −0.28, t = −3.14, P < 0.01; Fig. 1D; see Table 1). Between the two rank abundance parameters (mean and SD), we found that experimental warming marginally significantly increased the SD of the rank abundance log-normal distributions (β = 0.153, t = 1.61, P = 0.05), whereas no significant plant species richness effects were found (Table 1 and fig. S1). The mean of the rank abundance log-normal distributions was unaffected by both warming and plant species richness (Table 1 and fig. S1).

Fig. 1

Effects of experimental warming on nematode communities across the gradient of plant species richness.

(A) Nematode taxa richness. (B) Effective number of species. (C) Pielou’s evenness. (D) Rarefied taxa richness. Amb, ambient. Please see Table 1 for the details of the results.

Table 1

Results for nematode communities’ responses to experimental warming and plant species richness based on linear mixed-effects models.

Statistical significance is based on Wald type II χ² tests. Significant effects (P < 0.05) are given in bold. All regression coefficients are based on rescaled response variables. Please find the details of nematode responses in table S4 (mean responses with SD and number of samples). PSR, plant species richness; ENS, effective number of species; MTD, mean taxonomic distance.

Nematode responses	PSR			Warming			PSR × warming
	β	t	P	β	t	P	β	t	P
Richness	−0.012	−0.917	0.153	−0.158	−1.749	0.965	0.024	2.477	0.013
ENS	−0.010	−0.831	0.020	−0.233	−2.752	0.486	0.030	3.244	0.001
Pielou’s evenness	0.001	0.076	0.247	−0.200	−2.167	0.036	0.009	0.962	0.335
Rarefied richness	−0.003	−0.238	0.120	−0.282	−3.141	0.004	0.015	1.591	0.111
MTD	0.001	0.097	0.630	−0.137	−1.441	0.011	−0.005	−0.526	0.598
Mean (rank abundance distribution)	−0.004	−0.315	0.486	0.057	0.607	0.534	−0.002	−0.232	0.816
SD (rank abundance distribution)	−0.005	−0.416	0.317	0.153	1.612	0.052	−0.003	−0.319	0.750

The standardized effect size of mean taxonomic distance (MTD), which is calculated from the taxonomic distance matrix (see Materials and Methods for details), between taxa within plots decreased from a mean value slightly above 0 at ambient temperature to values less than 0 in warmed plots independent of plant species richness (β = −0.13, t = −1.44, P = 0.01; Fig. 2A). This indicates that experimental warming acted as an environmental filter for nematode communities across the gradient of plant species richness. Nonmetric multidimensional scaling (NMDS) analyses further revealed a marginally significant separation in nematode community composition (based on abundances) between the ambient and warmed treatments (F1,39 = 1.69, P = 0.05; Fig. 2B). We did not find a significant difference in nematode community composition between plant monocultures and 16-plant species mixtures (F1,39 = 1.23, P = 0.25).

Fig. 2

Effects of experimental warming and plant species richness on nematode taxonomic redundancy and community composition.

(A) Standardized effect sizes for the MTD of nematode communities in response to experimental warming across the gradient of plant species richness. (B) NMDS ordinations illustrating differences in nematode community composition among ambient and the highest warmed (~+3°C) plant monocultures (Mono) and ambient and the highest warmed (~+3°C) 16-plant species mixtures (Mix).

To test the generality of our results, we also explored the effects of other plant diversity indices (effective number of plant species and phylogenetic plant diversity) on nematodes in separate models. The effects of those two metrics of plant diversity were similar to those of plant species richness on nematode communities. For instance, neither effective number of plant species nor phylogenetic plant diversity showed significant effects on the MTD of nematode species; but we found similar interactive effects between either of these two plant diversity metrics and experimental warming on the effective number of nematode species as was found for plant species richness (table S2).

Nematode density decreased with warming in plant monocultures, whereas warming increased nematode density in the mixed plant communities, although this pattern was not significant [generalized linear mixed-effects model (GLMM); β = 0.01, z = 1.74, P = 0.08]. Among the four feeding groups (bacterial feeders, fungal feeders, plant feeders, and predators), we found a significant decline of predator density in warmed plots (~63%) independent of plant species richness (GLMM; β = −0.40, z = −1.58, P = 0.01). The density of the other feeding groups was not significantly affected by warming and plant species richness (table S3). Prey density (the sum of detritivore and herbivore densities) did not change significantly with warming and plant species richness (table S3). Further, we found a significant increase in prey richness in warmed and diverse plant communities, which was similar to the pattern for total taxa richness of nematodes (table S3). Predator richness did not change with warming and plant species richness (table S3).

Plant species richness and warming effects on plant biomass

Plant shoot biomass increased nearly threefold from plant monocultures to communities with 16 species (β = 0.03, t = 2.59, P < 0.001, log-transformed). For plant shoot biomass in the year of the nematode sampling (in 2013), we found no significant interaction effect between plant species richness and experimental warming (β = 0.008, t = 1.87, P = 0.06, log-transformed). Root biomass also substantially increased (nearly fourfold from plant monocultures to communities with 16 species; see table S4) with plant species richness (β = 0.05, t = 5.49, P < 0.001, log-transformed) in 2012 (no root data were available for 2013). However, root biomass was not affected by experimental warming (β = 0.10, t = 2.08, P = 0.09, log-transformed). We did not find any significant interaction effect of plant species richness and experimental warming on root biomass (β = −0.007, t = −1.28, P = 0.19, log-transformed).

Associations of plant biomass and predator density with nematode diversity

Predator density [potentially reflecting top-down effects; (, )] was positively associated with the effective number of nematode species in plant monocultures (Table 2 and Fig. 3A). The same pattern was found between predator density and MTD, but only in 16-plant species mixtures (Table 2 and Fig. 3D). Moreover, we found a significant negative interaction effect of plant shoot biomass [potentially reflecting bottom-up effects; (, , )] and predator density on the MTD of nematode communities in plant monocultures (Table 2). Root biomass was negatively associated with the effective number of nematode species in 16-plant species mixtures (Table 2). Most of the other associations between indicators of top-down and bottom-up effects and the two nematode diversity responses were statistically nonsignificant (Table 2).

Table 2

Regression results for the associations of bottom-up and top-down effects with the two nematode diversity measurements (ENS and MTD) using mixed-effects models in plant monocultures and mixed plant communities.

Statistical significance is based on Wald type II χ² tests. Significant associations (P < 0.05) are given in bold. All regression coefficients are based on rescaled response variables. S, shoot biomass; R, root biomass 2012; PD, predator density.

	Plant monocultures						Mixed plant communities (16 species)
	ENS			MTD			ENS			MTD
	β	t	P	β	t	P	β	t	P	β	t	P
Indicators of bottom-up effects
S	−0.252	−1.542	0.123	−0.083	−0.493	0.621	0.306	1.498	0.134	0.089	0.413	0.679
R	−0.091	−0.541	0.588	−0.172	−0.925	0.354	−0.485	−2.570	0.010	−0.017	−0.081	0.935
Indicator of top-down effects
PD	0.418	2.722	0.006	0.157	0.943	0.345	−0.05	−0.250	0.802	0.539	2.939	0.003
Interactive effects
S*PD	0.505	0.894	0.371	−1.331	−2.221	0.026	−0.934	−1.653	0.098	0.543	1.066	0.286
R*PD	0.126	0.427	0.669	−0.029	−0.09	0.928	0.327	0.327	0.528	0.653	1.345	0.178

Fig. 3

Relationships between predator density and nematode diversity measures in plant monocultures (A and B) and 16-plant species mixtures (C and D).

Statistical details of these relationships are provided in Table 2. Only significant (P < 0.05) relationships are shown by the dashed lines. These observed relationships do not imply causality (that is, predator density was not experimentally manipulated) but highlight the need to study the roles of trophic interactions in determining biodiversity in future global change experiments.

DISCUSSION

Our study has two key findings. First, warming increased nematode diversity (taxa richness and ENS) in the complex soil environment, as represented by high plant species richness, whereas it decreased nematode diversity in the simpler soil environment of plant monocultures. Second, warming appears to have acted as an environmental filter by selecting for more closely related (taxonomically) nematode species in both simpler and complex soil environments. Together, these results suggest that warmer environments can alter the community structure of nematodes by increasing the number of species () and their taxonomic redundancy in complex environments. Further, our results indicate that warming-induced shifts in top-down controls within the nematode community (that is, reduced predator density) and plant species richness–induced bottom-up controls (that is, greater plant biomass) associate with the observed nematode community responses in warmer environments.

The observed decline in nematode taxa richness and diversity in warmer plant monocultures could be due to decreased resource availability, lower soil water content, and homogeneous soil environment, as we hypothesized (, ). Warmer plant monocultures showed a slight increase in plant biomass production, but this effect was much lower than that of warmer species-rich plant communities (). This could mean that warmer plant monocultures were unable to subsidize warming-induced increases in resource demands of nematode communities. However, rarefied taxa richness and species evenness of nematodes declined because of warming alone (that is, independent of plant species richness). Therefore, our results not only indicate the inability of plant monocultures to sustain a greater diversity of nematodes in warmed conditions but also suggest warming-induced compositional shifts of nematode communities independent of plant species richness. Warming-induced responses of nematode communities may be driven not only by the available resources but also by processes such as physiological adjustments by species (, ).

Previous studies have also reported compositional shifts in nematode communities due to warming, that is, the abundance of some species decreased, whereas that of other species increased in response to warming, leading to altered nematode community structure (, , ). Negative effects of warming on nematode communities are often attributed to environmental stress, such as drying of the soil (, ). Water stress is also likely to become more pronounced in plant monocultures where a sparse vegetation (due to lower plant cover) is unable to reduce evaporation from the soil surface (, ). Our soil moisture data show greater soil water content in species-rich plant communities than in plant monocultures (fig. S3), whereas long-term data from this experiment show that soil water content decreases consistently with experimental warming independent of plant species richness (). Predatory nematodes are likely to be more sensitive to variations in soil moisture than prey nematodes because of their larger body size (table S3), which requires a thicker water film/larger water-filled soil pores for movement and feeding than is needed by smaller individuals. Studies have also shown that the community composition of nematodes varies with plant identity and composition, indicating that both plant-related resource diversity and availability are crucial for determining nematode community structure (, ). Hence, the variability in plant biomass production between plant monocultures (slight increase) and species-rich plant communities (stronger increase) due to warming () led us to expect that these variations in bottom-up forces may associate with the observed nematode diversity patterns in warmer environments. Although our analyses indicated the importance of bottom-up effects in driving nematode community structure, warming-induced shifts in predator density (potentially reflecting top-down effect) were also associated with nematode diversity patterns (Table 2 and Fig. 3).

Top-down forces help maintain greater species richness via the regulation of abundant and competitively superior prey species (, , ). The observed positive association between the predator density and the effective number of nematode taxa in plant monocultures provides some support for this theory (Fig. 3 and Table 2). A decrease in predation pressure in plant monocultures due to a drier and resource-limited environment (lower plant production) might have further magnified competitive hierarchies within the nematode communities, leading to a greater loss of species. However, we did not find any significant association between predator density and nematode diversity in species-rich plant communities (Table 2). Hence, our results regarding a potential link between predator density and nematode diversity in species-rich plant communities are inconclusive. The inconsistent associations between predator density and nematode species diversity between plant monocultures and species-rich plant communities point toward context-dependent effects of predation on species diversity () and merit further explorations.

The warming-induced increase in taxonomic relatedness across the gradient of plant species richness indicates that warming acted as an environmental filter for the nematode taxa that share similar responses to warmer environments. We also found some indication of nematode community convergence in the NMDS ordination space in warmed 16-plant species mixtures (Fig. 2B). Environmental stress has been shown to increase the similarity in species composition by eliminating nontolerant species (). Greater taxonomic similarity among nematode taxa in warmer plots accordingly suggests a warming-induced selection for species with greater tolerance to warmer and drier soil conditions. It is further noteworthy that predator density was significantly positively associated with the standardized effect sizes for the MTD in mixed plant communities (Table 2). This indicates that greater top-down control may promote overdispersion within the nematode communities, which likely reflects competitive exclusion of species of greater similarity (). The presence of generalist predators has been proposed as a driver of overdispersion within ecological communities (). We speculate that the observed decline in predator density may have reinforced environmental filtering processes in communities exposed to warmer environments, particularly when reduced predator densities may not reinforce the exclusion of taxonomically similar prey species.

In conclusion, warmer and complex environments may support consumer communities with a greater species diversity. However, these consumer communities may encompass taxonomically redundant species independent of environmental complexity, potentially due to warming-induced shifts in top-down and bottom-up forces. Moreover, a consistent negative effect of warming on nematode diversity in plant monocultures indicates their greater vulnerability to environmental stress in resource-limited environments. Intense land use practices in managed ecosystems (for example, agricultural monocropping) often result in communities with a very low number of plant species (–). Our results highlight the need for maintaining complex environments, such as plant species-rich communities, to offset detrimental warming effects on consumer diversity, but also point to the limitations of these approaches in maintaining taxonomically distinct consumer communities in a warmer world.

MATERIALS AND METHODS

Experimental design

We conducted the warming and plant diversity experiment at the long-term ecological research station at Cedar Creek Ecosystem Science Reserve in Minnesota, USA. Warming experimental units are nested within the long-term ongoing plant diversity experiment at Cedar Creek, which has been running since 1994 (, ). The starting conditions of the long-term biodiversity experiment (also known as “BigBio”) were fairly uniform across plots, that is, the whole area is on deep sands (formerly an abandoned agricultural field), and the top ~6 cm of the soil was removed to further homogenize the soil before the start of the experiment in 1994. The combined plant species richness and warming experiment [named “Biodiversity and Climate” (BAC); (, )], which began in 2008 (with first measurements taken in 2009), comprises 14 plant monocultures, nine 4-plant species mixtures, and 9 16-plant species mixtures. Plant mixtures are composed of randomly selected species from the pool of 18 species of the main biodiversity experiment. The pool of plant species represents four different plant functional groups: C3 grasses, C4 grasses, legumes, and nonleguminous species (, ).

Three warming treatments were fully crossed with plant species richness plots: ambient, ambient ~+1.5°C, and ambient ~+3.0°C (, , ). The heating was carried out using infrared heaters, which have been running since 2008 (). Low warming plots were heated using 600-W heaters, whereas the high warming plots were heated using 1200-W heaters. Metal frames were used to suspend infrared heaters at a height above 1.8 m from the ground. The ambient temperature plots also consisted of the metal frame and shade but without heaters to control for any shading effect (). The realized temperature in the highest warming treatment (+3°C) was, on average, about 2°C greater than in the ambient plots (when averaged over 25 cm aboveground to 30 cm belowground); however, this increase depended on the growing season. The difference in the realized temperature between the ambient and the highest warming treatments progressively decreased from the beginning of the growing season in May to peak biomass in July (from about 2° to 1°C) [for details, see the study of Cowles et al. ()].

Nematodes

Nematodes were sampled in August 2013 during peak plant biomass. Six soil cores (20 cm deep and 3 cm in diameter) were taken from each plot and carefully pooled for nematode extraction. A fraction of the pooled soil (20 g fresh weight) was used to extract nematodes using a modified Baermann method (). Nematodes were collected after 48 hours of extraction, preserved in 4% formaldehyde, and counted using a dissecting microscope. About 100 individuals from each sample were identified at random using an inverted Leica DMI 4000B light microscope. Nematode identification followed the work by Bongers () and was carried out to genus level (adults and most juveniles) and family level (remaining juveniles) (see table S1).

The number of nematode taxa was used to describe the taxa richness per plot. Nematode diversity was estimated as the effective number of species (ENS) () as an exponential of Shannon-Wiener index [H′ = −Σp × ln (p), where p is the proportion of the ith taxon; ()]. ENS provides diversity estimations least affected by the rarity or commonness of species within communities (). ENS and taxa richness are equivalent when all taxa within a plot have the same number of individuals (). Here, ENS is the effective number of nematode taxa. Pielou’s evenness (J′) was calculated as H′/ln(S), where S is the number of taxa (). To account for the differences in nematode abundance among plots, we also rarefied nematode taxa richness using the rarefy function in the “vegan” package for R statistical software (). Further, on the basis of the genus-level information, nematodes were assigned to four feeding groups: predators, which feed on both detritivore and herbivore nematodes; bacterial- and fungal-feeding nematodes as detritivores; and root-feeding nematodes as herbivores (). Omnivorous nematodes were not found in our samples (table S1).

On the basis of the topology of the Linnaean taxonomic information available for nematodes (class, order, family, and genus; see table S1), we estimated taxonomic relatedness (that is, MTD) between nematodes within each plot. Comprehensive phylogenetic information including reliable branch length estimates is, to date, not available for hyperdiverse taxonomic groups such as nematodes. However, diversity estimates based on hierarchical taxonomic information are often strongly related to those calculated from dated phylogenies (, ). We calculated the MTD on the basis of the interspecies cophenetic distance matrix constructed from a dendrogram depicting the taxonomic hierarchy [R package vegan; ()].

To assess whether nematode taxa within communities were taxonomically more or less related than expected or simply reflect differences in the species richness among plots, we compared the observed MTD values to those obtained from 999 randomly generated communities. Random communities were generated by shuffling the tips on taxonomic distance matrix [R package “picante”; ()]. This null model maintains the (i) levels of species richness within plots and (ii) species’ occurrence frequency among plots. The standardized effect size was calculated according to Gotelli and Rohde (): MTD.ses = (MTD.obs/MTD.exp)/sd(MTD.exp), where MTD.obs is the observed MTD, and the MTD.exp and sd(MTD.exp) are the mean and SD, respectively, across the 999 random communities. MTD.ses values <0 indicate that nematode taxa within communities are taxonomically more closely related, that is, taxonomically less distinct, than expected, given the levels of species richness—a pattern referred to as taxonomic clustering. MTD.ses values >0 indicate that taxa within communities are less related, that is, taxonomically more distinct, than expected—a pattern referred to as taxonomic overdispersion. Taxonomic clustering is often interpreted as a signature of (i) filtering processes that select for closely related species that share similar adaptations to a common environment () or (ii) competitive hierarchies that select for closely related, competitively superior taxa (). Taxonomic overdispersion could indicate competitive exclusion of closely related species ().

We also estimated aboveground biomass of plants in August 2013 using clip strips of 0.1 m × 1 m size from each plot. Roots were obtained in August 2012 (1 year before the nematode sampling) from soil cores of 5 cm (diameter) × 30 cm (deep) and later gently rinsed on a 1.5-mm mesh screen to remove soil. Plant shoot biomass and root biomass were determined after drying for >72 hours. Based on the species-specific plant biomass [from the study of Cowles et al. ()], we calculated the effective number of plant species (exponential of Shannon-Wiener index; see above) and used this as an alternative plant diversity measure to test whether other plant diversity metrics also yield similar results to that of plant species richness. We also calculated phylogenetic diversity of plant species across the richness gradient as an additional metric of plant diversity. Plant phylogeny was obtained from Zanne et al. (), and the phylogenetic diversity metric was based on the work of Faith () and calculated using the picante R package (). Hence, we were able to compare the effects of three metrics of plant diversity—plant species richness, effective number of plant species, and phylogenetic diversity of plants—on nematode communities in warmer environments.

Data analysis

The effects of experimental warming and plant species richness on nematode taxa richness, diversity, evenness, predator and prey richness, ratio between predator and prey density, and taxonomic relatedness were analyzed using linear mixed-effects models (Gaussian errors) with plots (warming plots nested within plant species richness plots) as a random intercept effect. Estimates of regression coefficients were obtained using the R package “lme4” (). Count data, such as nematode density (nematode counts per 20 g of fresh soil) and density of feeding groups, were analyzed using GLMMs with negative binomial error. Warming and plant species richness effects on plant biomass (shoot biomass and root biomass, separately) were analyzed using linear mixed-effects models with Gaussian errors with plots as a random intercept effect. We separately used effective number of plant species and phylogenetic plant diversity as two other plant diversity metrics as a covariate, together with warming treatments using the same mixed-effects models. These two additional metrics of plant diversity were used to explain the variations in two key nematode diversity metrics: effective number of nematode species and taxonomic relatedness of nematode species. We used Wald type II χ² tests to calculate the P values from the mixed-effects models using the “car” package (). Visual inspection of residual versus predicted values (that is, no correlation between the residuals and fitted parameters of the model) confirmed that the linearity assumptions were met for all the linear and GLMMs ().

Further, to compare the patterns of nematode community composition across warming and plant species richness, we compared the rank abundance distributions of nematode communities for each plot. Nematode abundance data were fitted to the Poisson log-normal distribution using the rad.lognormal function in the vegan package for R statistical software (). The two fitted parameters from the Poisson log-normal distribution (μ and σ) were then tested for warming and plant species richness effects using mixed-effects models. μ and σ were only estimated for plots with at least three taxa present because the plots with less than three taxa were insufficient for model fitting. We also used NMDS, based on Bray-Curtis dissimilarities, to visualize differences in nematode community composition. We used only plant monocultures and 16-plant species mixtures for NMDS analyses (R package vegan) in the lowest and highest warming treatments to visualize the extremes of warming and plant species richness effects. Further, we ran permutational multivariate analysis of variance (PERMANOVA) to test whether nematode community composition differences were significant using the adonis function in the vegan package.

Finally, we used regression models (linear mixed-effects, plots as random intercept) to associate interactive effects of potential bottom-up and top-down forces on nematode diversity patterns. Plant shoot biomass and root biomass were used separately as the bottom-up forces, whereas predator density was used as the top-down force. Plant biomass often correlates with microbial biomass in the soil (), which is the main basal resource of the soil food web. Predator density is a crucial top-down force to structure their prey communities (, ). Although predator diversity is a commonly used metric to relate predation effects on species diversity (, ), we chose predator density effects because warming-induced changes in predator densities were more pronounced than the changes in predator diversity (see results of nematode responses). Predatory nematodes can exert strong top-down control on prey nematodes (such as microbial-feeding nematodes) with important implications for prey community structure as well as cascading effects on microbial communities ().

Informed by the results of the regression models for plant species richness and warming effects on nematode communities (as explained above), we chose two nematode diversity measures for relating them with bottom-up and top-down effects: ENS and the metric for MTD. We ran these regression models for plant monocultures and 16-plant species mixtures separately to present a comparison for bottom-up and top-down effects on nematode diversity in simpler versus complex soil environments. All linear model assumptions were met with these models. Regression coefficients were estimated using the lme4 package (). All continuous variables for linear-mixed effects model (only those with Gaussian error terms) in this study were rescaled to improve the interpretability of regression coefficients () using the “rescale” function in the arm package ().