Nematodes in a polar desert reveal the relative role of biotic interactions in the coexistence of soil animals

Abiotic factors are major determinants of soil animal distributions and their dominant role is pronounced in extreme ecosystems, with biotic interactions seemingly playing a minor role. We modelled co-occurrence and distribution of the three nematode species that dominate the soil food web of the McMurdo Dry Valleys (Antarctica). Abiotic factors, other biotic groups, and autocorrelation all contributed to structuring nematode species distributions. However, after removing their effects, we found that the presence of the most abundant nematode species greatly, and negatively, affected the probability of detecting one of the other two species. We observed similar patterns in relative abundances for two out of three pairs of species. Harsh abiotic conditions alone are insufficient to explain contemporary nematode distributions whereas the role of negative biotic interactions has been largely underestimated in soil. The future challenge is to understand how the effects of global change on biotic interactions will alter species coexistence.


Supplementary Figure 2
Statistical evidence for negative biotic interactions requires three bits of evidence (Solid, thicker blue arrows): 1, verification that species co-occur non-randomly and, specifically, that they segregate; 2, statistical control of all other biotic and abiotic factors; 3, modelling of spatial autocorrelation.  Tables   Supplementary Table 1 Binomial GLMM of the probability of occurrence of Eudorylaimus. Spatial autocorrelation in the residuals (Figure 2) was modeled using an exponential function.

Extraction of Fauna
Nematodes, tardigrades and rotifers were extracted from soils using the modified sugar-centrifugation technique of 1 . In this method, fauna in soil are mixed with water and poured through stacked 500 and 38µm sieves to remove large rocks and debris. The resulting slurry is collected in a clean tube, centrifuged to separate the animals and soil from the water, and the supernatant removed. The pellet of animals and soil are then remixed in a dense sugar solution that is again briefly centrifuged. The soil and debris form a pellet and the animals are contained in the supernatant. The supernatant is washed over a 26µm sieve and specimens are collected, and stored in water. The animals are moved to counting dishes and identified using an inverted compound microscope. An advantage of the method is that it allows the collection of both living and dead animals. Nematodes were identified to species level 2-5 . Tardigrades and Rotifers were identified no further than phylum.
Soil microarthropds in Antarctica live in the first few cm of soil and indeed mostly between rocks and the substrate. Thus, the best method to represent their fauna if collecting them from the underside of small, flat, rocks [6][7][8] . In order to have quantitative estimates of their densities, specimens were collected within a 20m radius of where the soil sample for mesofauna was taken. Identification of specimens following the basic descriptions of 9 and follow up publications: see the general review by 7 for a full list of taxonomic references.

Molecular Methods
Total environmental DNA was extracted from soil samples using a CTAB protocol 10 modified for the X-tractor Gene liquid handling robot (Corbett Life Sciences, Concorde, NSW, Australia). Briefly, 0.7 g of soil was added to a microcentrifuge tube containing 0. Another surrogate for microbial biomass, total soil ATP, was measured in duplicates (triplicates where the duplicates did not match) using a 3M Clean-Trace Beverage Test Kit (Acorn Scientific, Auckland, NZ) with a modified protocol. In short, 100 µL of Extractant Buffer was added to 100 mg of soil and allowed to incubate for 60 seconds. 75 µL of ATP Assay Solution was then added to the sample, which was immediately read using a 3M Clean-Trace NG Luminometer (Acorn Scientific). Total soil ATP levels were recorded as The total number of true peaks in each sample was taken as a measure of taxon richness.
Soil respiration was measured by incubating 20 g (dry weight equivalent) samples of soils at 10°C for 26-28 days in an miniaturized respiromotric chambers 14 and the CO2 determined periodically by gas chemoatogrpahy (Varion90 GC fitted with a thermal conductivity detector) as described by 15 .

Modelling approach and rationale
Our aim was to detect significant negative covariation between pairs of species potentially competing for resources. We had to test for non-random patterns in species cooccurrence and assess whether the negative correlation between species co-vary with other biotic and abiotic factors (Supplementary Figure 2). Models must also account for spatial autocorrelation, which can be caused by unmeasured population processes such as dispersal dynamics (e.g., source/sink dynamics and historical legacy) as well as unmeasured abiotic and biotic factors.

Models and analyses (additional information on GLMMs)
Binomial GLMMs were applied to the presence/absence distribution of species, which was modeled as a function of abiotic predictors such as topography and soil properties (moisture, pH) and biotic factors such as the biomass of other taxa. Given our a-priori hypothesis, the main aim of these models was to detect a negative effect of Plectus on Scottnema (and vice versa) and neutral effects of Eudorylaimus on Plectus and Scottnema but we needed to ensure that these negative effects were not indirectly caused by other abiotic and biotic covariates.
We used GLMMs 16 to model spatial autocorrelation in the residuals using an exponential function, which was chosen after preliminary analyses. Specifically, to investigate spatial autocorrelation, we first ran a standard GLM that had the same structure of the GLMMs except for the spatial autocorrelation component. Residuals from these models were mapped and tested for autocorrelation using Moran's I statistics, autocorrelograms, and relevant tests

Abiotic predictors
We split abiotic factors into topographical factors, which describe the geomorphology of the study area, and soil factors, which account for microscale variability in soil factors such as moisture, pH and salinity. Principal component analysis (PCA) of the correlation matrix of topographical variables (see Supplementary Figure 1, PCA_top) showed that at least three axes were necessary to account for more than 2/3 of variance (Supplementary  Figure 3, PCA top) accounts for negative covariation between slope and the distance to the coast. Thus, in our linear models (see Table 1  Also in the case of the soil data matrix (Figure 3, PCA.soil), three PCA axes were necessary to meet the criterion of 2/3 of total variance: the first axis (PC1.soil, 33 %) accounts for a gradient in organic C and total N; the second axis (PC2.soil, 18 %) accounts for negative covariation between wetness or water and pH, while the third axis (PC3.soil, 18 %, not shown) is a gradient in electrical conductivity. In our linear models (see Table 1 and 2 in main text) we thus define PC1.soil, PC2.soil and PC3.soil as Organic matter gradient, Moisture gradient, and Salinity gradient respectively.

Biotic predictors
PCA of biotic variables (Figure 3, PCA.bio) identified the following three axes: the first one (PC1.bio, 33 %) accounts for a gradient separating sites with high microbial richness and high animal and microbial biomass from sites with low microbial richness and low animal and microbial biomass; the second one (PC2.bio, 18 %) accounts for negative covariation between, on the one hand, microbial biomass plus tardigrades and rotifers biomass, and on the other hand microbial richness and arthropod biomass; the third one (PC3.bio, 12 %; not shown) accounts for the fact that there were sites for which high biomass of mites and collembolans negatively correlated with the other biotic variables. . In our linear models (see Table 1 and 2 in main text) we thus define PC1.bio, PC2.bio and PC3.bio as