Small herbivores and abiotic heterogeneity promote trait variation of a saltmarsh plant in local communities

Study site

A natural successional gradient is present in the back-barrier saltmarsh of the island of Schiermonnikoog (53°30′N, 6°10′E), the Netherlands. Because this island expands eastward, thus the eastern part of the island is younger relative to the western part (Olff et al., 1997). The western part of the saltmarsh is enclosed for cattle grazing, while the eastern part of it is grazed by wild small herbivores such as spring staging geese, year-round present hares and rabbits. Hares and geese are the most abundant herbivores (Van De Koppel et al., 1996; Van Der Wal, Kunst & Drent, 1998; Van Der Wal et al., 2000a; Kuijper & Bakker, 2005; Schrama et al., 2015). Rabbits are very rare. Kuijper & Bakker (2005) found that biomass removed by rabbits in 2000 was 6.33% and 0.16% of the total biomass removed by small herbivores at the early and intermediate successional stage, respectively.

We used long-term hare and goose exclosures that were initiated in 1994 (details in Chen et al., 2019b). We selected exclosures located at the early and intermediate successional stages, which are approximately 2.5 km apart (Fig. 1A). Vegetation succession has undergone around 30 and 60 years at the early and intermediate successional stage, respectively. Age of the vegetation succession was counted from the year of vegetation establishment to 2016. The year of vegetation establishment was determined by checking aerial photographs (Olff et al., 1997). Grazing pressure from hares and geese was ca. 1.5 times higher at the early successional stage than the intermediate successional stage based on year-round dropping count in 2000 and 2016 (see Table S1 for more details; Kuijper & Bakker, 2005; Chen et al., 2019b; Chen et al., 2019a).

Hare and goose exclosures (one per each successional stage) were located in a similar elevation (early stage: 1.42 ± 0.004; intermediate stage: 1.44 ± 0.004; mean ± 1se; N = 24; m+NAP; Normal Amsterdam Water Level). Exclosures (8 m × 12 m and 6 m × 8 m at the early and intermediate successional stages) were made by chicken mesh (mesh width 25 mm) supported by wooden posts every 3.5 m to exclude hare grazing inside the exclosures. Exclosures were around 1 m in height, ropes were suspended on top of the wooden posts to stop geese flying into the exclosures. At the beginning of the exclosure experiment (1995), E. atherica rarely occurred (<2.5%; percent cover) inside and outside the exclosures at these two successional stages. In 2016 (22 years later), vegetation composition differed in the grazed areas and inside the exclosures (Fig. 1; also see Table S2 for species composition for the three most abundant species). Specifically, Artemisia maritima was the most dominant plant in the grazed area at these two successional stages, while E. atherica was the most dominant plant inside the exclosures.

Experimental design

To make our results spatially comparable, we sampled within ca. 6 m × 8 m inside and outside the exclosures at both successional stages in June 2017. We chose this size because this is the largest area we can sample inside the exclosures at the intermediate successional stage. The distance between the area and the exclosure measured ca. 80 m. We randomly set up 7 plots (1 m × 1 m) inside these two grazed areas and exclosures (Fig. 1). The distance between these plots ranged from 0.5 m to 9 m. We divided each plot into 25 grids (0.2 m × 0.2 m). Within each grid, we collected one individual stem of E. atherica (with roots), usually in the middle of the grid (Fig. 1). Elytrigia atherica did not occur everywhere, particularly not in the grazed area, thus sample size per plot varied from 9 to 25.

Trait measurement

We measured traits for individual stems of E. atherica in the field. We sampled individual stems without obvious grazing marks. We cut individual stems at the ground level and measured maximum height (cm) from the base to the highest point. We measured the maximum leaf width and length of the first fully grown leaf, usually the third leaf from the top. We also recorded whether individual stems flowered or not. After these measurements, each stem was stored in a paper bag, sealed, labeled, and then put in a self-sealing plastic bag to reduce water loss in tissues. At the end of each day, samples were brought to the lab, we measured fresh biomass of the individual stems and fully-grown leaves. We also measured their dried biomass (g) after oven-dry (70 °C) to constant mass (ca. 3 days). We calculated specific leaf area (mm² mg⁻¹) as (leaf width × leaf length)/dried leaf biomass. We calculated leaf dry matter content (mg g⁻¹) as dried leaf biomass/fresh leaf biomass. Specific leaf area and leaf dry matter content were not measured in the standard way (Pérez-Harguindeguy et al., 2016), and caution should be taken when comparing our data with other studies. Due to that reason, we focus on height, aboveground biomass, and flowering in the main text, we present the results for specific leaf area and leaf dry matter content in the Supplemental Information.

Genotyping and genotype richness and diversity

We used oven-dried leaf samples (ca. 2 mg per individual stem) for DNA extraction. We first shredded leaf samples into smaller pieces using tissuelyser. We then extracted DNA from each leaf sample using the CTAB method (Doyle & Doyle, 1987) and stored DNA samples at −20 °C before PCR. We amplified DNA using PCR with fluorescence-labeled primers. We used five microsatellite markers (ECGA89, WMS6, WMS44, WMS2, and ECGA89) originally designed for the other Poaceae species, Elymus caninus (Sun, Salomon & Bothmer, 1998) and Triticum aestivum (Röder et al., 1998). These five markers have been used for genotyping E. atherica in this system (Bockelmann et al., 2003; Chen, 2020). The PCR products from primer ECGA89, WMS6, and WMS44 were pooled together, while the PCR products from WMS2 and ECGA89 were pooled together. Additionally, 1 µL rROX was added in each sample as the internal size standard (Gene ScanTM–350 ROXTM, Applied Biosystem). We visualized the pooled PCR products using the 3730 DNA analyzer and scored the microsatellite peak patterns (height > 100) manually using GeneMapper. In total, we successfully genotyped 579 individual stems of E. atherica, but 2 were excluded for further analyses as some trait data measured in the field for these samples were missing.

Abiotic variables

We measured clay thickness and topographic variation for each plot (3 replicates). We measured clay thickness using a 2 cm Ø  soil corer with tick marks as a proxy for soil total nitrogen (Olff et al., 1997). Previous studies show that clay thickness is strongly positively correlated with soil total nitrogen in this system (i.e., total nitrogen (g) = 186 + 43 clay thickness (cm); r = 0.9, p < 0.0001; (Olff et al., 1997)). Also, clay thickness is used as a proxy of soil fertility (e.g., Schrama et al., 2017). We measured topography (elevation) using Trimble R8 (precision for elevation ca. 1 cm).

Calculating trait means, ITV, topographic variation, and variation in clay thickness

Although our focus is on ITV, to compare whether the effects of small herbivores on ITV are stronger than trait means, we also looked at trait means. We calculated means for each trait by averaging trait values over 9 -25 individual stems in each plot. We calculated ITV in each trait each plot as the standard deviation /mean. Topographic variation and variation in clay thickness were calculated as standard deviation /mean (over 3 samples) of elevation and clay thickness, respectively.

Calculating genotypic richness, genotypic diversity, and genetic differentiation

To calculate genotype richness and diversity per plot, we first calculated pairwise genetic distance using Dice dissimilarity from the R package ade4 (Dray & Dufour, 2015) based on the presence/absence matrix of 42 allele bands from those five markers. We then assigned genotypes based on dice dissimilarity, using the function “assignClones” from the R package polysat (Lindsay, Clark & Clark, 2018). We calculated genotypic richness as the number of unique genotypes detected divided by the number of individual stems genotyped for each plot. Genotypic diversity—taking into account the abundance of different genotypes—was calculated using the function “genotypeDiversity” with the index of “Shannon” from the package polysat (Lindsay, Clark & Clark, 2018). Small herbivores may select for some particular genotypes that are more resistant or tolerant to herbivore grazing (Kotanen & Bergelson, 2000), which may also impact traits. Therefore, we explored genetic differentiation using principal coordinates analysis (PCoA) from the R package ade4 (Dray & Dufour, 2015). We looked at genetic differentiation at the treatment level, that is, grazed and ungrazed across plots at each successional stage, result can be found in Fig. S1.

Effects of small herbivores on trait means and ITV

We used analysis of variance, function “lm”, to look at the effects of small herbivores on means and ITV of each trait at the early and intermediate successional stage, separately. We checked model assumptions by visually inspecting residual plots for homogeneity and quantile–quantile plots for normality. In all cases, models fitted the data reasonably well.

Direct and indirect effects of small herbivores on ITV

We estimated the direct, indirect effects through altering genotypic richness and diversity, indirect effects through altering abiotic variables, and total effects (sum of direct and indirect effects) of small herbivores on ITV. We estimated these effects for each trait at each successional stage based on the standardized path coefficients using structural equation models from the R package “lavaan” (Rosseel, 2012). Note that leaf dry matter content was not evaluated in structural equation models due to limited sample size (see online supporting text). Grazing (0: ungrazed, 1: grazed), genotypic richness, genotypic diversity, topographic variation, and variation in clay thickness were included in the models. Note that genotypic diversity was not included for ITV at the early successional stage, as it was significantly correlated with genotypic richness (variance inflation factor > 4). See online supporting text for an example of a structural equation model for ITV in height at the intermediate successional stage and calculation of the direct and indirect effects. Effects are significant when p ≤ 0.05, while marginally significant when p > 0.05 and p < 0.1. Data analysis was performed in R3.5.3 (R Core Team, 2020).

Effects of small herbivores on trait means and ITV

At the early successional stage, small herbivores significantly decreased means in height, aboveground biomass, and flowering of E. atherica. Also, small herbivores significantly increased ITV in flowering and marginally increased ITV in height, but had no effects on ITV in aboveground biomass. At the intermediate successional stage, small herbivores significantly decreased means in height, aboveground biomass, and flowering. However, small herbivores had no effects on ITV in these traits (Fig. 2; see Table S3 for test statistics). Small herbivores did not impact mean and ITV in specific leaf area and leaf dry matter content at both successional stages except that small herbivores decreased mean in specific leaf area at the intermediate successional stage (Table S3).

Direct and indirect effects of small herbivores on ITV

The structural equation model reveals that overall (summing up the direct and indirect effects), small herbivores tended to promote ITV in height and flowering, but not in other traits at the early successional stage. While, small herbivores tended to promote ITV in height, this did not occur in other traits at the intermediate stage. Moreover, these positive overall effects of small herbivores on ITV in height and flowering were mainly attributable to their direct effects. Although the overall effects of small herbivores on ITV in aboveground biomass and specific leaf area were not significant, small herbivores also directly significantly impacted ITV in these two traits at the early successional stage (Fig. 3; Fig. S2). See Table 1 for the standardized path coefficient estimated from SEM for the direct, indirect, and total effects of small herbivores on ITV.

At the early successional stage, small herbivores did not have significant effects on genotypic richness and diversity, but genotypic richness decreased ITV in flowering. At the intermediate successional stage, small herbivores did not have significant effects on genotypic richness, but significantly reduced genotypic diversity. However, genotypic richness and diversity did not impact ITV in any traits at the intermediate successional stages (Fig. 3; Fig. S2).

Small herbivores did not impact topographic variation and variation in clay thickness at both successional stages. However, these variations significantly increased ITV in all traits at the early successional stage except that the effects of topographic variation on ITV in flowering were not significant. The effects of these variations on ITV were less apparent at the intermediate successional stage except that variation in clay thickness marginally decreased ITV in height (Fig. 3; Fig. S2).

Effects of small herbivores on trait means and ITV

Our results suggest that small herbivores could impact both means and variations in key functional traits of E. atherica, but their effects differed in different traits. At the early successional stage, we found that small herbivores significantly decreased means in three traits measured (height, aboveground biomass, and flowering) and they significantly promoted ITV in two traits (height and flowering). At the intermediate successional stage, small herbivores significantly decreased means in four traits (height, aboveground biomass, flowering, and specific leaf area) and they marginally promoted ITV in one trait (height; Fig. 2; Fig. S2). This suggests that trait means and variations may be driven by different processes. Moreover, we only observed ITV in height and flowering may be because these two traits have higher plasticity than other traits, and traits with a higher level of plasticity show higher ITV (Givnish, 2002). Furthermore, small herbivores had similar effects on trait means at both successional stages, but had stronger impacts on ITV at the early successional stage. This may be because the effects of small herbivores on trait means and variations were mediated by community composition in the field. At the early successional stage, the preferred plant, F. rubra, was more abundant in the grazed area (ca. 3 times relative to that of the intermediate successional stage), thus small herbivores (although more abundant) may graze more on F. rubra rather than E. atherica. At the intermediate successional stage, the abundance of F. rubra deceased, while the abundance of E. atherica increased in the grazed area (ca. 5 times relative to that of the early successional stage), thus small herbivores (less abundant) may also substantially grazed on less preferred E. atherica. Therefore, overall, the effects of small herbivores on trait means were similar at these two successional stages. Meanwhile, non-preferred species (A. maritima) was more abundant in the grazed area at the early successional stage (ca. 2 times relative to that of the intermediate successional stage; see Table S2 for detailed percent cover of these species). It may facilitate young E. atherica escape from herbivore grazing. Therefore, ITV was more apparent at the early successional stage than the intermediate one. Previous studies in this system have shown that hares and geese are the important drivers for plant community composition and structure along this successional gradient (Olff et al., 1997; Kuijper & Bakker, 2005; Chen et al., 2019b). Here, extending the current knowledge, we show that small herbivores may also impact trait means and variations, which in turn could impact community composition and structure (Whitlock, Grime & Burke, 2010).

Direct and indirect effects of small herbivores on ITV

As we hypothesized, the direct effects of small herbivores through selective grazing increased ITV. This is not only true for height and flowering but also for aboveground biomass and specific leaf area, despite the overall effects of small herbivores on these two traits were not significant. As explained in the previous paragraph, the direct effects of small herbivores may be mediated by community composition especially the proportion of preferred and non-preferred plants. Herz et al. (2017) found that local neighborhood diversity can explain a large amount of ITV in German meadows and pastures, possibly through increased plant-plant interactions (e.g., competition). In this study, plant diversity was higher in the grazed than the ungrazed plots at both successional stages after 22-year grazing (Table S2). Thus, higher plant diversity may also contribute to increased ITV in height under grazing at both successional stages. Future studies/experiments looking at the effects of herbivores on ITV in plant communities with the same composition are needed to fully separate the effects of selective grazing from the effects of neighboring plants.

We found that small herbivores only decreased genotypic diversity at the intermediate successional stage. A previous study in the west part of this saltmarsh found that cattle grazing does not impact genotype richness and diversity of this plant, but cattle grazing alters its morphological traits in the field such as decreasing height and leaf width (Veeneklaas et al., 2011). Here, we found that genotypic richness and diversity generally had no significant effects on ITV. Microsatellite markers are commonly used in studying plant genetics (Jarne & Lagoda, 1996; Vieira et al., 2016), but because they are selectively neutral, they may not control gene expression for the traits measured here. More generally speaking, the genetic basis of most traits in many wild plants is unknown. Thus, it is difficult to know in advance which genes or genomic regions to consider, which further hinders inferring a causal relationship between genetic variation and phenotypic variation (Merilä & Hendry, 2014). A low correlation between variation in morphological traits and DNA markers was also reported before (Kolliker et al., 1998). Additionally, we found genetic differentiation of E. atherica in the grazed and ungrazed areas at the early successional stage, but not at the intermediate successional stage (Fig. S1). However, more data (presumably from more exclosures) are needed to consolidate this result. Thus, the observed trait changes (at both successional stages) may not be associated with genetic changes in E. atherica.

In contrast, we found no effects of small herbivores on topographic variation and variation in clay thickness within plots. These abiotic variations may be primarily induced by flooding and inundation in this saltmarsh, which may overrule the effects of small herbivores. However, variations in these two abiotic variables had substantial positive impacts on ITV particularly at the early successional stage. Abiotic environment may be more heterogenous at the early successional stage because this stage has more frequent sediment deposition than the intermediate successional stage (Schrama, Berg & Olff, 2012). A previous study conducted in the western part of this saltmarsh also found that topographic variation at the small spatial scales (0.1 m² and10 m²) is positively correlated with species richness and cattle grazing additionally increased species richness (Ruifrok et al., 2014). Thus, abiotic variations, even at very small spatial scales, may play an important role in promoting ITV and altering other plant community properties.

Long-term exclosures

Although we used 22-year old herbivore exclosures, our data cannot answer whether duration of herbivore grazing impacts ITV, as we only collected data for one year. To our knowledge, no studies have compared the effects of short- and long-term grazing on ITV. Didiano et al. (2014) found that tolerance to rabbit grazing decreased as the age of the exclosures increased in F. rubra, the most abundant plant in Silwood Park, England. Therefore, grazing duration may also impact ITV.