Nematodes and the effect of seasonality in grassland habitats of South Africa

Abstract Nematodes in South Africa have mainly been studied for their diversity and agricultural importance. However, the ecological status of nematodes and the effect of seasonal variation in local grasslands remain unknown. For this reason, a nematode study was conducted in the Telperion Nature Reserve and represented the first ecological study in a natural grassland area in South Africa. In total, 104 soil samples were collected during four consecutive seasons from 2015 until 2016 in three habitats, viz. (i) open grassland, (ii) shrubland with rocky outcrops, and (iii) riparian zone. From these the nematode community structure and soil ecosystem status were studied. In total, 93 genera from 50 families were recorded with herbivores and bacterivores being the most abundant trophic groups in all three habitats. Linear mixed models revealed that season had an overwhelmingly dominant impact on the condition, food web status, and functioning of the soil ecosystems with pairwise comparisons indicating that significantly higher values were recorded during winter. Interestingly, this seasonal shift can largely be attributed to fluctuations in the populations of only a few nematode groups (namely Aporcelaimellus, Dorylaimidae, Iotonchus, and Mononchus) with high colonizer-persister values. Although the reason for the higher abundance of specific nematode groups recorded during the winter is not explicitly clear, it is possibly linked to reduced competition from other soil fauna. This study clearly shows that further investigations are required to better understand the dynamics of grassland ecosystems.

and erosion control (Hewins et al., 2018;Bengtsson et al., 2019;Gibson and Newman, 2019). However, although the structure and ecological status of grassland nematode communities have been investigated in some parts of the world (Ekschmitt et al., 2001;Biederman and Boutton, 2009;Kergunteuil et al., 2016), a literature survey revealed no ecological reports on the nematode communities of South African grasslands.
Studying the nematode communities of grasslands is not only relevant from an ecological perspective, but also from a conservation perspective. Nematodes are valuable indicators of soil ecosystem disturbance affected by, for example, agricultural activities and climate change (Ferris and Bongers, 2009;Zhong et al., 2017;Zhang et al., 2020). Therefore, studying the nematode communities of natural, undisturbed grasslands will generate valuable baseline data for monitoring disturbance and implementing conservation policies. This would be especially valuable for South Africa as the country's grasslands are under continuous threat from overgrazing, cultivation, and urban expansion (Haddad and Butler, 2018). Woody (or bush) encroachment, driven by rising temperatures and increasing atmospheric carbon dioxide levels, also presents a major threat (Skowno et al., 2017). This while only 2.0 to 2.8% of South Africa's grasslands are protected within conserved areas (Carbutt et al., 2011;Haddad and Butler, 2018).
When studying terrestrial ecosystems, it is also important to consider seasonal variation as water availability, temperature, and vegetative growth can have an important effect on the structure and functioning of soil ecosystems (Bongers and Ferris, 1999;Cardoso et al., 2016;Song et al., 2016;Siebert et al., 2020). Song et al. (2016), for example, showed that higher precipitation rates significantly increased nematode abundance and richness in a grassland ecosystem. Bakonyi et al. (2007), in turn, found that the effects of soil drying and warming on nematode communities were dependant on the presence of vegetation with the smallest effect recorded in grasses. Ultimately, understanding how seasonal variation affects soil ecosystems is necessary in order to accurately assess and monitor the potential threats posed by anthropogenic activities and climate change.
For these reasons, this study was undertaken and aimed at (i) studying the nematode community structure in the open grassland, shrubland with rocky outcrops and riparian zone habitats of the Telperion Nature Reserve and (ii) investigating the effect of seasonal-induced changes on the condition, food web status and functioning of the soil ecosystems.

Site description
The study was conducted in the Telperion Nature Reserve (Mpumalanga province, South Africa) ( Fig. 1) that has a total surface area of 9,061.3 ha. This reserve is located in the Rand Highveld Grassland (of which only 1% is conserved) and is dominated by grassy plains interspersed with rocky outcrops and woody plant species (Grobler, 1999). The greater region is characterized by strong summer rainfall (570-730 mm per annum) and dry and cold winters (frost may occur frequently). This is a fire prone ecosystem with high lightning flash densities making lightning-induced fires common. The lithology of Telperion Nature Reserve is dominated by Arenite-Conglomerate, producing dystrophic or mesotrophic soils with some red soils, as well as rocky areas with miscellaneous soils (Grobler, 1999). Three streams flow through the reserve that originates from higher lying wetlands and sponge areas.

Rainfall and temperature data
Climatic data (rainfall and temperature) from January 2014 to December 2016 were obtained from the South African Weather Service (www.weathersa. co.za). The total amount of rain for 90 days prior to each sampling event (as listed in the following section) was calculated, while the mean, minimum, and maximum daily temperatures were extracted for the same period. Since sampling was based on climatic factors (as listed in the following section) and thus not evenly dated over 1 year, this 90-day period was considered before each sampling event.

Site selection and sample collections
Within Telperion Nature Reserve three main terrestrial habitats, namely, (i) open grassland (OG), (ii) shrubland with rocky outcrops (SRO), and (iii) riparian zone (RZ) were identified for investigation. Sites within each habitat were selected based on accessibility as some areas were restricted and/or not accessible by vehicle or foot. The elevation of the sites in the each habitat ranged as follows: OG: 1333 to 1482 m, SRO: 1361 to 1480 m, and RZ: 1304 to 1383 m above sea level. Ultimately, 26 sampling sites ( Fig. 1) were selected and included: 11 OG, 10 SRO sites, and 5 RZ sites.
In total, 104 soil samples (one per site per season) were collected from the three terrestrial habitats during four consecutive seasons [Winter (11 June 2015), Spring (30 November 2015), Summer (24 February 2016), and Autumn (19 April 2016)]. These sampling dates were selected based on climatic patterns through consultation with the reserve manager. Winter samples were the first to be collected when frost was at its peak. Spring sampling followed the first rains, which due to a drought in South Africa came later than normally expected. Summer sampling was undertaken when grasses reached maturity and produced seeds, while Autumn samples were collected before the first frost occurred.
Multiple discrete soil samples (500 cm 3 each) for the analysis of nematode community structure were collected from each site up to a depth of 20 cm using a garden trowel, which was cleaned between sampling sites to prevent contamination. Samples were collected beneath plants (if present). For the analysis of selected soil properties (soil texture and total organic carbon), a sub-sample was collected from each site. The latter was performed only during the summer (February 2016) sampling interval as the selected soil properties were unlikely to change significantly during the sampling period. All the collected samples were labeled and transported in cool boxes to the Nematology Unit of the Agricultural Research Council -Plant Health and Protection (Roodeplaat, South Africa). Samples were stored at 10°C until further analysis.

Physical and biological properties of soil samples
The selected physical and biological properties were analyzed by Eco-Analytica (North-West University, South Africa) as follows: total organic carbon (C) content of the soil samples were determined using the loss-on-ignition method (Donkin, 1991) and the soil texture as described by Laker and Dupreez (1982).

Statistical analysis
Rarefaction curves were used to compare taxa richness between habitats. This was achieved by calculating the sample-based Mao Tau estimator (with inter-and extrapolation) using Estimate S 9.1 Software Package. Furthermore, differences in the abundance of nematode trophic groups between habitats were investigated by pooling samples from across seasons. Abundance values were log 10 (x + 1) transformed and visualized, while statistical significance between trophic groups and habitats was inferred using two-way analysis of variance (ANOVA) and Tukey's multiple comparisons tests. These analyses were performed using Graphpad Prism 6 Software Package.
Selected nematode-based indices (i.e. maturity index 2-5, enrichment index, structure index, channel index, basal index, and metabolic footprints) were used to quantify the condition, food web status, and functioning of soil ecosystems (Ferris et al., 2001;Ferris and Bongers, 2009;Ferris, 2010). These indices were calculated using the Nematode Indicator Joint Analyses (NINJA) online tool (Sieriebriennikov et al., 2014). The Graphpad Prism 6 Software Package was used to create a violin plot of the maturity index 2-5, a measure of soil ecosystem condition. Also, values of the structure and enrichment indices were used to plot a faunal analysis and characterize the food web status of the studied soil ecosystems (Ferris et al., 2001). Furthermore, linear mixed models with pairwise comparisons were used to determine if the independent variables, i.e. season (Winter, Spring, Summer, and Autumn) and habitat (OG, SRO, and RZ) significantly affected (singularly and interactively) the calculated nematode-based indices. Prior to this analysis the nematode-based indices were log transformed (where needed). All linear mixed models were performed using SPSS Statistics 25 Software Package.
Lastly, a redundancy analysis (RDA) was used to investigate the relationship between the response (nematode-based indices) and explanatory variables [(factors (seasons and habitats) and numeric variables (soil texture and total organic carbon)]. Predictor effects of the listed explanatory variables were calculated using a Monte-Carlo permutation test (Šmilauer and Lepš, 2014). These multivariate analyses were performed and illustrated on a biplot using Canoco 5 Software Package. Significance for all univariate and multivariate analyses was regarded at p < 0.05.

Environmental conditions
The total rainfall for 90 days prior to each sampling event is illustrated in Figure 2A. This shows a general increase over time with the lowest and highest rainfall recorded prior to the Winter and Autumn sampling events, respectively. However, it is important to note that South Africa experienced a drought in 2015. According to the South African Weather Service (www.weathersa.co.za), the Telperion Nature Reserve region received only 388 mm from January to December 2015. By contrast, the recorded annual rainfall in 2014 and 2016 were 794 and 748 mm, respectively.
Temperature (mean, minimum, and maximum) values are illustrated in Figure 2B with the lowest and highest mean temperatures recorded, as expected, prior to the Winter and Summer sampling events, respectively. The largest range between minimum (−2.6°C) and maximum (31.5°C) values were recorded before Winter sampling. Finally, the soil texture and organic carbon content of the collected and analyzed soil samples are provided as Supplementary material (Table S1).

Nematode community structure
Nematode communities associated with the studied habitats were investigated in terms of their taxonomic composition and richness, as well as their trophic abundance. In total, 93 genera, representing 50 families, were recorded in samples collected from Telperion Nature Reserve (Supplementary  material Tables S2-S6) and were representative of the five major nematode trophic groups, namely, herbivores (Table S2), bacterivores (Table S3), fungivores (Table S4), omnivores (Table S5), and predators (Table S6). The largest number of taxa were found in the shrubland with rocky outcrops (SRO) habitat (77 genera), followed by the riparian zone (RZ) (74 genera) and open grasslands (OG) (68 genera) habitats. In terms of the number of recorded genera per trophic group, bacterivores dominated with 30 genera recorded in both the SRO and RZ habitats, while 24 bacterivore genera were found in the OG habitat. This was followed by herbivores A comparison of taxa richness between habitats was made using rarefaction curves (Fig. 3). These curves illustrate the sample-based observed (solid line) and extrapolated (dotted line) Mau Tau richness, while 95% confidence intervals are indicated as shaded bands. When comparing richness at the 20 samples mark (total number of samples collected  at the RZ habitat), it is clear that RZ presented the greatest richness with 84 genera. This was followed by SRO and OG, which presented 73 and 65 genera, respectively. Nonetheless, it is important to note that these curves did not reach an asymptote (leveling), suggesting that additional sampling might be required to record the complete nematode community. To this end, extrapolation revealed that more than 50 samples per habitat are likely required at which point an estimated 98, 91, and 83 genera would be predicted to be collected at the RZ, SRO, and OG habitats, respectively.
Lastly, the nematode trophic abundances per habitat were considered. Figure 4 shows that herbivores and bacterivores were the most abundant trophic groups, followed by fungivores, omnivores, and predators. The multiple comparisons tests revealed that although no significant differences occurred between herbivores and bacterivores, both these trophic groups differed significantly from fungivores, omnivores, and predators at all three habitats. Also, significant differences between fungivore, omnivore, and predator abundances were recorded at all three habitats with the exception of no significant difference between fungivores and omnivores/predators at the RZ habitat. No significant differences in nematode trophic abundances were recorded between habitats.

Soil ecosystem condition and food web status
The maturity index 2-5 was used to investigate the soil ecosystem condition of the studied grassland habitats. A violin plot (Fig. 5) illustrates the median, quartiles, minimum, maximum, and data distribution (shape outline) of this index per habitat per season. Minimum and maximum values ranged from 2 (RZ in Spring) to 3.7 (SRO in Winter), respectively, while the highest and lowest median values for all three habitats were recorded during the Winter and Spring, respectively. Also worth noting is that the SRO habitat had the highest median values during all four seasons. This habitat was thus classified as having the most stable soil ecosystem. The results from the linear mixed models (Table 1) showed that habitat and season had a significant effect on the maturity index 2-5. Yet, pairwise comparisons revealed that these effects were only evident between seasons, which showed significantly higher maturity index 2-5 values in Winter compared to the rest of the seasons. Differences between Spring, Summer, and Autumn were not significant.
The status of the soil food webs were investigated using the nematode faunal analysis (Fig. 6). This analysis showed that all three habitats had the highest structure (complexity) in the   Winter with the SRO and RZ habitats classified as maturing and enriched with balanced (bacterial vs. fungal) decomposition pathways. The OG habitat was classified as mature and fertile characterized by fungal dominated decomposition. In Spring, however, nematode community composition in all three habitats represented degraded and depleted food webs. During Summer and Autumn, the communities represented either degraded and depleted or mature and fertile food webs. The faunal analysis also revealed the clustering of habitats (i.e. low variability between habitats) especially during Winter and Spring. The effect of season and habitat on the two indices used to determine the food web status, namely the structure index and enrichment index (Fig. 6), were further investigated using linear mixed models. This revealed that while both season and habitat had a significant effect on the enrichment index, only season had a significant effect on the structure index. Pairwise comparisons showed, as with the maturity index 2-5, that communities in Winter had significantly higher enrichment and structure compared to the other seasons, while the OG habitat had significantly lower enrichment compared to the SRO and RZ habitats. Finally, a redundancy analysis illustrated on a biplot (Fig. 7) was used to investigate the relationship between the explanatory and response variables. The explanatory variables accounted for 42% (p < 0.001) of the observed variation in the response variables with axes 1 and 2 representing 35.3% (p < 0.001) and 5.5% (p < 0.001), respectively. This analysis revealed that the three habitats could be differentiated based on the predominant soil fraction, i.e. the OG, RZ, and SRO habitats consisted of more sandy, silty, or clayey soils, respectively. Nonetheless, only clay presented a significant effect of 4.3% on the response variables. In addition, organic carbon levels were positively correlated with the RZ habitat and enrichment index, but had only a minor effect (< 4%) on the response variables. Furthermore, the redundancy analysis biplot showed that a positive correlation existed between the OG habitat and the channel and basal indices. The SRO habitat, in turn, presented a positive correlation to the maturity 2-5 and structure indices. Finally, the Monte-Carlo permutation test revealed that 30% of the variation in response variables can be attributed to seasonality.

Soil ecosystem functioning
Metabolic footprints (Table 2) were used to assess the soil ecosystem functioning, i.e. the magnitude of functions and services delivered by the soil ecosystems. Generally, relatively high herbivore footprints were recorded, while especially high herbivore footprints were observed in Winter and Spring at the SRO and RZ habitats, respectively. However, even though clear differences were recorded in the herbivore footprint between seasons and habitats, the linear mixed models (Table 1) revealed that none of these differences were significant. Of the remaining footprints, significant seasonal effects were only recorded for the predator and structure footprints. Pairwise comparisons revealed that the predator footprints were significantly higher in the Winter compared to Spring and Summer, while the structure footprints were significantly higher in the Winter compared to Spring and Autumn. No significant habitat effects were recorded for any of the footprints.

Discussion
The nematode community structure Grasslands typically support a high richness of nematode taxa (Song et al., 2017). Yeates and Lee (1997) and Bell et al. (2005) reported 32 and 70 genera, respectively, from tussock grasslands in New Zeeland. Popovici and Ciobanu (2000), in turn, listed between 33 and 67 genera from different  grassland regions in Romania, while Čerevková (2006) reported 64 genera in grasslands from the Slovakia Republic. The present study, however, reported a substantially higher number of nematode genera of 93 in total. Rarefaction curves suggest that increased sampling efforts at these sites would result in further taxa being added to that total. At the RZ habitat, for example, extrapolation of the Mau Tao richness indicator suggests that close to a 100 nematode genera could be recovered with increased sampling efforts. The high number of recorded and estimated genera is important as it is commonly accepted that a high richness in nematode taxa are indicative of healthy and functioning soil ecosystems (Ekschmitt et al., 2001;Ferris, 2010;Sánchez-Moreno and Ferris, 2018). A high richness of herbivore taxa in natural environments are often the result of a high diversity in plant species (Hodda et al., 2009;Cortois et al., 2017). According to Hodda et al. (2009), some evidence suggests that a potential reason for this is the lack of competition resulting from root systems creating very defined niches.
Another important consideration when looking at nematode community structure is the abundance distribution between trophic groups (Ferris et al., 2001;Ferris, 2010;Van Den Hoogen et al., 2019). A recent study by Van Den Hoogen et al. (2019) investigated the abundance of nematode trophic groups across all major terrestrial biomes and showed that, on average, bacterivores dominated, followed by herbivores, fungivores, omnivores, and predators. The same trend was observed for temperate grasslands, which is the classification under which the South African grassland biome is listed (Van Den Hoogen et al., 2019). However, even though bacterivores are reportedly the most abundant trophic group on a global and biome scale, studies have reported regional grassland nematode communities with large and even dominating herbivore populations (Porazinska et al., 2003;Ayres et al., 2008). Similarly, during the present study, no significant difference was recorded between the two most abundant groups, namely herbivores and bacterivores, at any of the habitats.
It is worth noting that herbivores of economic importance (e.g. Meloidogyne) occurred in relatively high numbers at some sampling sites. Other herbivores (e.g. Helicotylenchus, Rotylenchulus, and Scutellonema) were present in all the habitats during all four seasons. The presence of these herbivore taxa, sometimes in high numbers, poses a threat to crop production in this region. According to various studies conducted in production areas where major grain crops (e.g. maize and soybean) are produced, Meloidogyne is generally the most predominant nematode pest. Substantial yield losses in maize (up to 60%) (Riekert and Henshaw, 1998) and soybean (up to 100%) (Smit and De Beer, 1998;Fourie et al., 2013) due to infection by this genus demonstrate its adverse impact on crop production. Moreover, Helicotylenchus, Rotylenchulus, and Scutellonema are also abundant in local crop fields, indicating their potential impact on crop production (Fourie et al., 2001;Bekker et al., 2016;Mbatyoti et al., 2018). Therefore, grasslands may act as reservoirs for herbivores that pose a threat to crop production should such areas be converted into agricultural fields.

Seasonal and habitat induced changes in grassland ecosystems
Seasonality had a large and significant effect on the condition, food web status and functioning of the studied grassland soil ecosystems. This effect was most evident progressing from Winter to Spring. During Winter the habitats were enriched with stable ecosystems and structured food webs, while depletion with reduced stability and structure were evident in Spring. Furthermore, during the Winter increased ecosystem functioning were evident with significantly higher predator and structure footprints. This even though low temperatures and reduced precipitation (as recorded during the 90 days period prior to Winter sampling) typically exert a negative influence on nematode communities (Verschoor et al., 2001;Song et al., 2017;Andriuzzi et al., 2020). According to Song et al. (2017), the optimal temperature for the survival and reproduction of soil nematodes ranges from 20 to 25°C, while extremes below 5 and above 30°C significantly inhibit development. Similarly, water availability influences primary production and thus energy flow into the soil food web, while also regulating organic carbon decomposition (Andriuzzi et al., 2020). It should be noted that the especially high herbivore footprints can be attributed to the large number of Meloidogyne that were present in some of the samples. Metabolic footprints are calculated using the biomass of adult females (Ferris, 2010), which can substantially increase footprint values if large numbers of, for example, some sedentary endoparasitic nematodes are recorded.
In contrast to seasonal effects, habitats and the measured soil properties had minimal influence on the studied soil ecosystems. This despite habitat type, soil texture and organic carbon levels being known to substantially influence soil faunal communities (Du Preez et  Nonetheless, the positive correlation (as evidenced by the redundancy analysis) between organic carbon and the enrichment index is in line with the general understanding that higher organic carbon levels support opportunistic bacterial-feeding nematodes (Sánchez-Moreno and Ferris, 2018). Increased levels of soil organic carbon at the RZ habitat is therefore likely the reason for this habitat also having the highest nematode community enrichment levels during all four seasons in the soil food web analysis. Conversely, the basal and channel indices were positively correlated with the OG habitat, indicating that soil ecosystems associated with this habitat were likely more degraded, but with greater fungal decomposition (Ferris et al., 2001;Ferris, 2010;Sánchez-Moreno and Ferris, 2018).
Interestingly, the abundance of nematode taxa (Supplementary material: Tables S2-S6) made it clear that the dramatic shift in nematode community in Winter can largely be attributed to fluctuations in the populations of only a few predacious nematodes (namely Aporcelaimellus, Iotonchus, and Mononchus), as well as one omnivorous nematode group, the Dorylaimidae. These taxa are assigned a colonizer-persister (c-p) value of 4 (Dorylaimidae and Iotonchus) and 5 (Aporcelaimellus and Mononchus) (Sieriebriennikov et al., 2014), which means that they are regarded as persisters (or k-strategists). Their presence in large numbers therefore infers greater ecosystem stability and food web structure (Ferris et al., 2001;Sánchez-Moreno and Ferris, 2018), as was recorded during Winter at the SRO and OG habitats. The reason for the increased number of members of these nematode groups recorded during Winter is not explicitly clear, however, findings from previous studies provide some potential answers. According to De Ley (1992), desert and dune sands often contain a remarkably high proportion of dorylaims, while Arpin (1969) repeatedly found dorylaims (predators and omnivores) in desiccation experiments. De Ley and Mundo-Ocampo (2004), in turn, stated that the ecological sensitivity of these nematodes may have been overestimated and that certain species may be capable of surviving desiccation and freezing. Omnivorous nematodes typically also have versatile feeding habits and can probably interact at various levels of the soil food web (Hanel, 2003). It is therefore possible that dorylaims in the studied Rand Highveld Grasslands of South Africa are not adversely affected by the conditions present in Winter and may even increase in numbers following reduced competition for resources from other faunal groups. The possibility that certain nematode families and genera are favored by Winter conditions (low temperature and rainfall), infers that nematode community responses to altered conditions are taxa dependent (see also Papatheodorou et al., 2004).

Conclusion
The grassland habitats of Telperion Nature Reserve host a diversity of nematode taxa, many of which are likely still to be observed or described. Although the studied nematode communities were dominated by herbivores and bacterivores, it was the less abundant nematode groups, i.e. omnivores and predators, which had a significant influence on the condition, food web status, and functioning of the soil ecosystems. Furthermore, evidence suggests that only a few taxa are largely responsible for seasonal shifts in measurable soil ecosystem parameters.