Experimental treatments
WENNDEx manipulates three global change factors (nighttime warming, increased winter precipitation, and N addition) in a fully crossed, completely randomized 2 X 2 X 2 design. WENNDEx has five replicates of each treatment combination for a total of 40 plots of 3.0 × 3.5 m each. All replicate plots contain both blue and black grama grass. Nighttime warming was imposed using lightweight aluminum fabric blankets drawn across each warmed plot at night, ~0.5 m from the ground to trap outgoing longwave radiation, which increased winter nighttime air temperature by an average of 1.1 °C and summer nighttime air temperature by 1.5 °C (Collins et al. 2017). Dataloggers controlling shelter movements retracted the blankets when wind speeds exceeded a threshold (to prevent damage) and when rain or snow occurred. Based on long-term climate records, El Niño rains increase average winter precipitation in our area by 50%; more frequent and intense El Niño events are predicted by climate models (Collins et al. 2017; Cai et al. 2015a,b). From 2006 2019 we supplemented winter precipitation each year using an irrigation system with reverse osmosis (RO) water because groundwater did not match rain chemistry. Rain was added in six experimental events each winter (January-March) to mimic the total amount (50 mm) and size distribution (four 5 mm events, one 10 mm event, and one 20 mm event each winter) of typical winter-storm events during El Niño years. To augment nitrogen, we used watering cans to add 2 g N m-2 y-1 as NH4NO3 prior to the monsoon season because NH4-N (57%) and NO3-N (43%) contribute approximately equally to N deposition at SNWR (Báez et al. 2007). Control plots received the same amount of RO water, the equivalent of a 2 mm rain event. In addition to the imposed treatments, on 4 August 2009, a lightning-initiated fire began on the Sevilleta National Wildlife Refuge. By 5 August 2009, the fire had reached WENNDEx, which was burned extensively though not entirely. Approximately 50% of plots burned, and those plots which did not burn were burned within three weeks by the US Fish and Wildlife Service fire crew. Prescribed fire is separately used as a grassland management tool at the Sevilleta National Wildlife Refuge with >10 year return interval. All plots were burned, so it is unlikely that fire affected our results and would not have affected treatments differentially.
Soil collection and processing
We collected soil nematodes at two time points to capture potential seasonal differences in their sensitivity to environmental changes. On each sampling date, we took two soil cores (2.5 cm diameter, 20 cm depth) from each of the 40 WENNDEx plots, near a grass-dominated area toward the center of each plot. We combined the two cores for data collection. Samples were collected during the spring (26 June 2020) and during fall (10 September 2020). Each core was collected in the early morning, coincident with sunrise, and placed into a labeled plastic bag for transportation to the laboratory. Samples were processed within 10 h of collection by placing a subsample from each bag into a modified Baermann funnel (Ingham 1994). Briefly, we placed a thin layer of 50 g of soil on tissue paper that acted as a filter mounted above water in wire-mesh pans (2 mm). We allowed the water to contact the soil via filter paper on the wire mesh. Taking advantage of the hydrophobicity of nematodes, nematodes drilled through the filter paper into the water. After 48 h of incubation at 22 °C, the nematodes in the water were concentrated on a 500-mesh sieve (25 μm aperture). After the total number of nematodes was counted, 100 specimens per sample were randomly selected and identified to genus. If there were fewer than 100 specimens per sample, all specimens were identified to genus (Yeates et. al 1993). Morphological identification was completed under a microscope at 400X magnification (eyepiece × objective, 10 × 40).
Statistical analyses
Hypothesis 1: Environmental changes interact to alter nematode diversity, abundance, community composition, and functional groups more than single factors alone. We calculated diversity metrics (genus richness, Shannon diversity index (H’), and the inverse of Simpson’s evenness) from the matrix of counts of identified nematode genera from each plot and sampling date using the vegan package in R (Oksanen et al. 2013). We also summed nematode genera in each of the following functional groups representing trophic guilds: bacterivores, fungivores, herbivores, or omnivore-carnivores (Yeates et. al 1993). For nematode diversity metrics, total abundance and functional group abundance, we built general linear mixed effects models with all three treatments, their interactions, the repeated effect of time (2 sampling dates) and the random, repeated effect of plot nested with warming × precipitation × nitrogen. Mixed models were constructed using lmer in the lme4 package in R (Bates et al. 2015) and met assumptions of homogeneity of variances and normality of residuals. We present analysis of deviance results of likelihood ratio tests using the Anova function in the car package (Fox and Weisberg 2018). We decomposed treatment interactions using post-hoc Tukey HSD tests among pairs with the emmeans package and provided corresponding P-values in the results using pairwise comparisons (Lenth 2021).
To detect treatment effects on nematode community composition, we used the matrix of counts of identified nematode genera to test for warming, nitrogen, and water treatment effects and all possible treatment interactions using a repeated measures design with the random effect of plot nested in the warming × precipitation × nitrogen treatment via perMANOVA in Primer 6.0 (Clarke and Gorley 2006). Models used 100 random restarts and 9,999 permutations. We examined whether dispersion differed among treatment combinations using permDISP, and we identified which genera contributed most to differences among treatments with SIMPER analyses using 9,999 permutations, also in Primer 6.0 0 (Clarke and Gorley 2006).
Hypothesis 2: Nematode community responses to environmental change track the abundance of plant biomass or plant community composition. First, we used a general linear model to correlate total nematode abundance with total estimated plant biomass from Baur et al. (2021), including the fixed factor of sampling date and the date × plant biomass interaction to test for seasonality in the relationship. Second, we examined relationships between the matrix of nematode genus-level composition and the matrix of plant species composition (via estimated biomass) both for the season of collection and for the season prior to nematode collection, in order to detect lagged effects of plants on nematodes. These analyses used Mantel tests for matrix correlations, implemented in Primer (Clarke and Gorley 2006). Black grama grass (Bouteloua eriopoda) is the foundational plant species in this ecosystem, we used DIST-LM to associate genus-level nematode composition with black grama biomass and linear models to associate the abundance of each nematode functional group with black grama biomass, including biomass during the season of collection and also during the prior season in order to detect potential lag effects, using Primer for DIST-LM and lm in R for regression analysis (Clarke and Gorley 2006). We also conducted similar DIST-LM for summed biomass of annual forbs, perennial forbs, and all grasses. In prior results, the largest plant responders to interactive treatments were the annual forbs (Collins et al. 2017).
