Sexual reproduction trait expressions of grassland species along a gradient of nitrogen: phosphorus stoichiometry

Plant investment in sexual reproduction is affected by nitrogen (N): phosphorus (P) stoichiometry. It has been suggested that an important adaptation to strong P limitation is reduced investment in sexual reproduction. We aim to investigate the specific influence of N:P on sexual reproduction performance within and between grassland species. Eleven grassland species were selected in ten plots covering N limitation, co-limitation and P limitation. Nutrients in soil and above-ground biomass were determined, plus soil pH and soil moisture. A range of sexual reproduction traits were measured as a proxy for investment in sexual reproduction. At the intraspecific level: compared with N-limited plots, in P-limited/co-limited plots, flowering time was later, flowering period in individuals was shorter, and number of flowers (inflorescences) per individual was smaller. At the interspecific level, in P-limited/co-limited plots, species had a significantly earlier flowering time and a longer seed stalk and seed panicle, than those in N-limited plots. However, flowering period was shorter and number of flowers (inflorescences) per individual was smaller under P limitation/co-limitation. Moreover, significant correlations between soil pH and soil moisture, and sexual reproduction performance of the selected grassland species were also found. P limitation/co-limitation restrict the sexual reproduction of grassland species, which may hamper their dispersal capacity. We recommend future studies further analyze the relationship between soil pH and N:P stoichiometry and the influence of soil pH, as well as soil moisture on sexual reproduction performance of grassland species in addition to analyzing N:P stoichiometry.

groundwater (Dentener et al. 2006;Vitousek and Howarth 1991). P enrichment, however, may even be more disruptive than that of N, due to increased mining of P-containing rocks for fertilizer production (Falkowski et al. 2000). Nowadays, after decades of heavy application, most soils in Europe are saturated with P (Obersteiner et al. 2013). N and P enrichment has been changing not only the absolute availability, but also the relative availability of N and P globally. The latter one is reported with dramatic impact on various terrestrial ecosystems, and started to draw increasing attention recently (Wassen et al. 2005;Fisher et al. 2012;Fujita et al. 2014).
Compared to other ecosystems, European grasslands have a rich flora and may develop a high species diversity (Pärtel et al. 1996;Schaminée et al. 1996). Many of these grasslands are naturally N-limited and their native flora has evolved mechanisms to cope with low N availability (Chapin 1980). It is therefore not surprising that the prevailing paradigm is that N enrichment diminishes or removes N limitation and enables fast-growing grasses with high N requirements to outcompete slow-growing forb species, which leads to loss of diversity (Bakker and Berendse 1999;Bobbink et al. 2010;Stevens et al. 2010). However, the type of nutrient limitation has rarely been analyzed and empirical data from a wide array of European grasslands have revealed that N limitation, P limitation or N and P co-limitation occur frequently, whereas K limitation is very rare (Wassen et al. 2005). Wassen et al. (2021) found that, in line with stoichiometric niche theory, these stoichiometric ratios influence species pools at the continental level and especially affect the distribution of threatened species.
These effects of stoichiometric ratios have been attributed to adaptations of plant species to environments of low nutrient availability or to specific types of nutrient limitation, and to trade-offs with growth and competitive ability (Fujita et al. 2014;Grime 2001;Roeling et al. 2018). Several plant traits have been identified that enhance nutrient acquisition, improve nutrient retention and/or nutrient use efficiency, and together contribute to a conservative nutrient strategy (Aerts and Chapin, 2000;Grime 2001;Lambers et al. 2008). These traits have been used to explain the ability of species to cope with low N availability (symbiotic N 2 fixation, chitinase production: Aerts and Chapin 2000;Leake and Read 1990); or low P availability (organic acid and phosphatase production and exudation: Fujita et al. 2014;Lambers et al. 2008;Wang and Lambers 2020), but irrespective of the type of nutrient limitation, most traits (high root: shoot ratio, formation of cluster roots, resorption from senescing tissues) may be beneficial (Aerts and Chapin 2000;Lambers et al. 2008). A competition experiment for two grasses along an N:P stoichiometric gradient could only partly explain the competitive advantages of plant traits (Olde Venterink and Güsewell 2010). Under changing environmental conditions, a species' plasticity in its traits becomes a key determinant for its future performance (Berg and Ellers 2010;Callaway et al. 2003), which leads us to suggest that phenotypic plasticity for acquiring and using nutrients is beneficial in a world in which nutrient availabilities are in flux. For one particular trait (phosphatase production), Fujita et al. (2010) showed that N:P stoichiometry is the most important determinant for such plastic responses. Moreover, Fujita et al. (2014) showed that plants adapted to P-limited environments invest less in sexual reproduction. Specifically, by comparison with common species, species of P-limited environments have lower seed number and seed investment, start flowering later, have a shorter flowering period, and have a longer lifespan (perennials rather than annuals). The combination of these traits suggests that species adapted to P-limited environments by investing less in sexual reproduction. Furthermore, Wang et al. (2019) showed that grassland species can also respond to changes in absolute and relative P availability by exhibiting intraspecific plastic responses in investments in sexual reproduction, i.e. species invested less in sexual reproduction under P-limited conditions than they do under N-limited or co-limited conditions. Such adaptations may enable species to survive in low P environments, provided that (1) other more competitive species are unable to grow under these circumstances and (2) eutrophication by P does not occur. In the present study we focused on sexual reproduction traits and investigated whether they vary not only between species but also within species along an N:P gradient in the field.
We hypothesized that along this gradient: 1) Species' responses in terms of investment in sexual reproduction are plastic, i.e., under P-limited conditions they invest less in sexual reproduction than they do under N-limited conditions (intraspecific differences); 2) Characteristic species of P-limited environments invest less in sexual reproduction than do common species or species characteristic of N-limited conditions (interspecific differences). To test these hypotheses we analyzed the expression of five types of sexual reproduction traits of 11 selected grassland species along a gradient of relative N and P availability in natural herbaceous plant communities. We selected unfertilized grassland sites in nature reserves that taken together we expected to encompass an N:P gradient (De Mars and Garritsen 1997;Roeling et al. unpublished data;Wassen 2017). Traits of flowering phenology, flowering stalk and panicle length, percentage of flowered individuals, and number of flowers (inflorescences) per individual were selected as the proxies of sexual reproduction of the selected grassland species. In detail, first flowering date (FFD), peak of flowering (PF) and last flowering date (LFD) were included because it has been shown that nutrient fertilization in nutrient-poor conditions may advance flowering (Putterill et al. 2004). The duration of the flowering period in individuals (FPI) and in the population (FPP) were included because it has been shown that duration of flowering is related to investment in sexual reproduction (Fujita et al. 2014). Stalk height with seed panicle (SHP) and stalk length with seed panicle (SLP) were included because we assumed that the greater the height at which the seeds are released, the greater the chance that the seeds will be dispersed over longer distances (Soons et al. 2004). Seed panicle length (SPL) was included because longer seed panicles might indicate more seed production. Percentage of flowered individuals (PFI) was included to indicate the floral promotive effect along nutrient conditions (Krekule and Kohli 1981). Lastly, number of flowers or inflorescences per individual (NF(I)I) was included because in other studies (e.g., Burton 1943), traits like this appeared to be a good indicator of sexual reproduction.

Study area and plot selection
The study was carried out in moist non-fertilized grasslands in three nature reserves located in two regions c.42 km apart: the Middenduin (MD) nature reserve in the west of the Netherlands, near Haarlem (52°23ˊN 4°35ˊE), and Laegieskamp (LK) and Koeiemeent (KM) nature reserves, which are 40 km further east, near Naarden-Bussum (52°16ˊN 5°8ˊE) (Fig. 1). Syntaxonomically the studied grasslands belong to Molinio-Arrhenatheretea / Parvocaricetea (Westhoff and Den Held, 1969). The three sites were chosen because, taken together, they encompass a gradient from N to P limitation in which nutrient availability overall is relatively low. The two regions are subjected to similar N input from atmospheric deposition (1500-2000 mol/ha) (https:// www. atlas natuu rlijk kapit aal. nl/ groot schal ige-stiks tofde posit ie-neder land). All three nature reserves are mown annually in summer by a nature management organization, and the hay is removed. The two regions differ slightly in their major soil type: peaty sand (De Mars and Garritsen 1997) prevails in Naarden-Bussum, and sand (van Raalte 2014) in Middenduin. The altitude of both regions is between 0 and 1 m above mean sea level. Their mean annual precipitation is 887 mm and their mean annual temperature is 11.2 ℃ (https:// www. clo. nl/ indic atoren/ nl0004-meteo rolog ische-gegev ens-in-neder land? ond= 20883).
We selected a total of 10 2 × 2 m 2 plots ( Fig. 1): two in LK, three in KM, and five in MD. The location of each plot was recorded using GPS. Plots were chosen in homogeneous vegetation.

Field surveys and chemical analyses
Species composition, vegetative productivity, and soil sampling Vegetation relevées were made on 19 and 20 June 2017, at the peak of the growing season. All herbaceous species were identified on the basis of Meijden (2005), and the cover was estimated with the Braun-Blanquet (1964) scale (Table S1).
To obtain a proxy for site productivity, aboveground standing biomass was harvested on 19 and 20 June 2017, which were the peaks of the growing seasons of the year (Roeling et al. 2018). The total above-ground biomass of each plot was estimated by clipping the above-ground vegetation in four 20 cm × 20 cm squares, one at each of the four corners of each plot. Occasional shrubs and litter were excluded from biomass sampling. Samples were dried immediately at 75 °C for 48 h in the lab and weighed. Our estimate of productivity is an underestimate of annual net above-ground productivity since we did not include bryophyte production (cf. Pałczynski and Stepa 1991;Wassen et al. 1995).
In May 2017, 4 soil samples (upper 10 cm) were collected in the four corners of each plot, using a 2 cm diameter auger. The soil samples were collected in plastic bags and stored at 4 ˚C for a maximum of 48 h before analysis.

Determination of plant N, plant P and plant K
Vegetation material was ground to pass through a 0.5-mm sieve. Total plant N of dried plant material was measured with a C/N elemental analyzer (NA1500, Carlo Erba-Thermo Fisher Scientific); total plant P was measured by trace element analyzer (S2, PICOFOX, Bruker). Total plant K was determined spectrophotometrically with a Sherwood Model 420 Flame Photometer.

Determination of available soil N and soil P
To estimate available soil N, c. 4 g of moist soil per sample was placed into a 50 ml centrifuge tube (polypropylene) and extracted with 40 ml 2.0 M KCl by shaking the suspension immediately for 30 min at 160 rpm (Barrett et al. 2002). Next, samples were centrifuged at 4000 rpm for 4 min and then filtered through a 0.45 um nylon filter (National 25 ml Nonsterile Centrifugal Filters, Thermo Scientific). The filtrate was stored at -20 ℃ until further analysis. NH 4 -N and NO 3 -N concentrations in the extract were measured on a colorimetrical auto-analyzer (AA3, Bran & Luebbe). Soil N was calculated as mg N/g dry wt. soil, using the specific moist soil weight of each sample and soil moisture concentration.
To estimate available soil P, c. 0.4 g of moist soil per sample was placed into a 50 ml centrifuge tube (polypropylene) and extracted with 40 ml ammonium lactate-acetic acid by shaking the suspension for 4 h at 100 rpm (Houba et al. 1979). Samples were then centrifuged at 4000 rpm for 4 min and filtered through a 0.45 um nylon filter (National 25 ml Nonsterile Centrifugal Filters, Thermo Scientific). The filtrate was stored at 4 ℃ until further analysis. P concentration in the extract was determined with ICP-OES (Spectro Acros, Spetro). Soil P was calculated as mg P/g dry wt. soil, using the specific moist soil weight of each sample and soil moisture concentration.

Determination of other environmental variables
The investigation was completed by determining soil pH and soil moisture content, which are known to be relevant for grassland species richness and community composition (Cornwell and Grubb 2003;Duprè et al. 2010;Silvertown 1980).
• Soil moisture: c. 5 g of moist soil was oven dried at 70 °C for 48 h to determine gravimetric water content. • Soil pH: soil and water were mixed in a ratio 1: 2.5 (c. 4 g air-dry soil and 10 g water) and then shaken for 2 h. Soil pH was measured in suspensions (Houba et al. 1979).

Nutrient limitation definitions
We use the plant nutrient concentrations-referred to in this paper as plant N, plant P, and plant K-as indicators of available nutrient concentrations (Fujita et al. 2014). When plant N:P < 13.5, the corresponding plot was defined as N-limited; when 13.5 ≤ plant N:P ≤ 16, the corresponding plot was defined as N-and P-co-limited; when plant N:P > 16, the corresponding plot was defined as P-limited. We made this decision based on previous research showing that above-ground plant N:P ratio is a reliable indicator to infer the type of nutrient limitation of the corresponding plot (Güsewell and Koerselman 2002;Olde Venterink et al. 2003;Wassen et al. 2005). For K limitation and K and N co-limitation, we referred to the standard of Olde Venterink et al. (2003), i.e., N:K > 2.1 and K:P < 3.4. This method of measuring nutrient in the above-ground vegetative biomass is an integrative way with the advantage of measuring over the growing season compared to soil nutrient concentration, which only provides one snapshot over time (Fujita et al. 2014). Soil nutrient concentration in the current research was used as a complementary environmental factor apart from plant nutrient concentration.

Selected species and their flowering status
We monitored the flowering phenology and a set of traits representative for the investment in sexual reproduction of a preselected set of 11 species with relatively high occurrence in the 10 plots. All individuals of these species if occurring in a plot were monitored throughout the growing season of 2017. The species list and information about their flowering status are presented in Appendix 1 Table 4. In detail, in terms of high occurrence of flowered individuals per plot (more than 5), species Succisa pratensis, Ranunculus flammula, Anthoxanthum odoratum, Juncus acutiflorus, and Carex panicea were mainly in P-limited/co-limited plots, whereas species Lysimachia vulgaris, Plantago lanceolata, Holcus lanatus, and Parnassia palustris were mainly in N-limited plots. A high occurrence of flowered individuals per plot of the remaining two species, i.e. Lythrum salicaria, Lotus uliginosus, were found on both P-limited/co-limited plots and N-limited plots (Appendix 1 Table 4).

Investigation of sexual reproduction performance
The sexual reproduction performance of the selected species was monitored almost weekly from the middle of April to the end of August 2017. After the growing season, the frequency of investigation was gradually reduced from once/week to once/two weeks, and to once/month, ending at the end of November 2017. The variables and traits recorded in each plot are discussed below.

Number of budding, flowering, and faded individuals
We define a budding individual as an individual of a species with only buds at the time of recording, a flowering individual was an individual of a species with open flowers or inflorescences at the time of recording, and a faded individual was an individual of a species with only faded flowers or inflorescences at the time of recording. These three variables of each selected species occurring in each plot were counted as accurately as possible at every field visit, from the first flowering data (FFD) to the last flowering date (LFD) of that species (for definitions of FFD and LFD, see below). Using these three variables of each species, we generated an overview of the flowering phenology of that species over the whole growing season.

First flowering date (FFD)
We define the day on which the first flower or inflorescence of a certain species occurred in a certain plot as the first flowering date of that species in that plot.
Since the interval between consecutive visits was c. 1 week, if a species did not have flowering individuals in one visit but did the next time, the first flowering date (FFD) of that species in that plot was taken to be the median date of those two visits' dates. The date of the earliest flowering individual of all monitored species in all plots combined was assigned the value 0 (which was that of Anthoxanthum odoratum in plot KM3). For statistical analysis, we then converted the other FFDs, which were originally in Julian dates, into numbers (i.e., days after day 0).

Last flowering date (LFD)
We define the day on which the last flower or inflorescence faded in a certain plot as the last flowering date of that species in that plot. If necessary, as explained above, we used the median date of two consecutive visits. The LFD was converted into a day number relative to the earliest FFD, as mentioned above.

Peak of flowering (PF)
We define the day on which the maximum number of flowers or inflorescences of a certain species occurred in a certain plot as the peak of flowering of that species in that plot. To obtain this value, the number of flowers or inflorescences of each species occurring in each plot were counted as accurately as possible at every field visit from FFD until LFD.

Flowering period in individuals (FPI) and in the population (FPP)
In most studies, flowering duration is measured monthly or in seasons (Fujita et al. 2014;Lahti et al. 1991). However, in our previous greenhouse experiment (Wang et al. 2019), we observed that flowering duration measured in days provides much better data at the intraspecific level under various nutrient treatments. Given the lack of a commonly agreed or specific definition of flowering period, we used two indices to assess flowering period in days. 1) Flowering period in individuals (FPI): this method consisted of marking random flowers or inflorescences of a species per plot (how many flowers or inflorescences depended on the total number of flowers or inflorescences of that species in the plot and ranged between 1 and 8). A waterproof marker denoting the FFD of a certain flower or inflorescence was tied to the stalk of the flower or inflorescence on its first flowering day and the day on which that flower or inflorescence started to fade was recorded. The duration between these two dates was the FPI of that flower or inflorescence. The average of all the monitored flowers or inflorescences of a certain species in a certain plot was considered as the flowering period in individuals (FPI) of that species in that plot. However, the result of this trait measure was probably not accurate, since the flowering period of individual flower or inflorescence was short (mostly a few days) and there was an interval of c. one week between two consecutive visits.
2) Flowering period in the population (FPP): this method consisted of calculating the time in days between the first flowering date (FFD) and the last flowering date (LFD) of a certain species in a certain plot.

Stalk height with seed panicle (SHP) and stalk length with seed panicle (SLP)
We measured the stalk height with seed panicle relative to the soil surface (SHP). Next, after straightening the stalks manually so they were vertical, we measured the stalk length (SLP) by measuring the distance between the soil surface to the top of the seed panicle (Heady 1957). The first method is indicative of seed release height, while the second one (which is a more standardized measurement) is indicative of a plant's investment in the structure of seed stalks.

Seed panicle length (SPL)
Seed panicle length was measured from the bottom to the top of each seed panicle.

Percentage of flowered individuals (PFI)
We counted the number of budding individuals, flowering individuals, faded individuals, and vegetative individuals without sexual reproduction of each species in each plot in each time visiting. Percentage of flowered individuals (PFI) was calculated by dividing the largest number of the total of budding individuals, flowering individuals, and faded individuals, by the largest number of the total of budding individuals, flowering individuals, faded individuals, and vegetative individuals.

Number of flowers or inflorescences per individual (NF(I)I)
The total number of flowers or inflorescences (NF(I)) of a certain species was acquired by counting the number of all flowers or inflorescences of that species in a certain plot accurately during each visit. If a species had started to set seed, the seed panicles were counted as well. The largest number recorded was taken to be the NF(I) of that species in that plot. Number of flowers or inflorescences per individual (NF(I)I) was calculated by dividing NF(I) by the total number of flowered individuals, i.e., the total of budding individuals, flowering individuals, and faded individuals of each species in each plot.

Data analysis
One-tailed t-tests were used for comparing two data groups (P limitation/co-limitation and N limitation) with normal distributions. The variables investigated were plant variables (plant N, plant P, plant N:P, and plant biomass), and soil variables (soil N, soil P, soil N:P, soil pH, and soil moisture); and intraspecific and interspecific differences in trait expression between P-limited/co-limited plots and N-limited plots. We used one-way ANOVA for comparing three data groups (plots of site LK, plots of site KM, and plots of site MD), which included plant variables (plant N, plant P, plant N:P, and plant biomass), and soil variables (soil N, soil P, soil N:P, soil pH, and soil moisture). Games-Howell procedure was applied for post-hoc testing since homogeneities of variances were significantly different. Simple linear regression analyses were used to ascertain the relationships between soil nutrient concentrations (i.e., soil N, soil P, and soil N:P), and plant nutrient concentrations (i.e. plant N, plant P, and plant N:P) and plant biomass. Simple linear regression analyses were also used to identify relationships between plant N:P, soil N:P, soil pH, soil moisture, and first weighted then averaged traits, i.e., FFD, PF, LFD, FPP, FPI, SLP, SPL, PFI, and NF(I)I of all species in each plot. All the analyses mentioned above were performed in SPSS 23.0 (SPSS, Chicago, U.S.A.); the corresponding figures were also created in SPSS. A detrended correspondence analysis (DCA) was carried out using CANOCO 5, to explore the relationship between environmental variables and species composition; the corresponding figure was also created in CANOCO 5.

Nutrient conditions of the grasslands
The selected plots covered all three types of nutrient limitation (Fig. 2). Plant N:P ratios indicate that the Laegieskamp plots (LK1 and LK2) are limited by P (N:P > 16) and KM1 and KM2 are co-limited by N and P (N:P < 13.5). The other six plots are limited by N (KM3, MD1-5) (13.5 ≤ N:P ≤ 16) (Fig. 2). Examining the plot results after the plots had been divided into two nutrient limitation classes (P limitation/colimitation (N:P ≥ 13.5; i.e. LK1, LK2, KM1, KM2) and N limitation (N:P < 13.5; i.e. KM3, MD1, MD2, MD3, MD4, MD5)) revealed that the plots with P limitation/co-limitation had significantly higher mean plant N:P than plots with N limitation (Table 1). No K limitation or K and N co-limitation was found in the selected plots (data not shown).
Plant N were similar for all plots, with largest variation occurring in the Middenduin area. Although plant N were significantly higher in Laegieskamp and Middenduin than in Koeiemeent, the differences were small (in LK and MD, + 11.00% and 7.29% higher than KM; Appendix 2 Table 5). On the other hand, plant P of Laegieskamp was significantly lower than that in Koeiemeent and Middenduin, and the Critical thresholds for P and N limitation cf. Olde Venterink et al. (2003) are indicated with dotted horizontal lines. P limitation relative to N when N:P > 16, N and P co-limitation when 13.5 ≤ N:P ≤ 16, N limitation relative to P when N:P < 13.5 Table 1 Differences in plant N and plant P (mg/g), plant N:P ratio, above-ground living phanerogam biomass (g/m2), soil N and soil P (mg/g), soil N:P ratio, soil pH, and soil moisture (%) (mean ± SD) between P-limited/co-limited plots (N:P ≥ 13.5; i.e. LK1, LK2, KM1, KM2) and N-limited plots (N:P < 13.5; i.e. KM3, MD1, MD2, MD3, MD4, MD5) Significant differences are indicated by * (one-tailed t-test P < 0.01) P-limited/co-limited plots N-limited plots differences were large (in LK, 53.69% lower than in KM and 59.76% lower than in MD; Appendix 2 Table 5). In Koeiemeent, the above-ground biomass is significantly lower than in the other areas (Appendix 2 Table 5) and the P-limited/co-limited plots have lower biomass than N-limited plots (Table 1).
Soil pH and soil moisture of the grasslands The three areas differ significantly in soil pH, with Laegieskamp and Koeiemeent having acidic environments and Middenduin having alkaline conditions (Appendix 2 Table 5) . P-limited/co-limited plots had a significantly lower pH than N-limited plots (Table 1). Soil moisture was higher in Laegieskamp than in Koeiemeent and Middenduin throughout the season and was higher in P-limited/co-limited plots than in N-limited plots (Table 1, Appendix 2 Table 5, Appendix 3 Fig. 4).
Intraspecific and interspecific sexual reproduction performance of the selected grassland species along the gradient in N:P ratio Because of the absence of overlap in species composition between plots, only a limited number of the 11 species that were monitored (Table 2) were found with significant intraspecific differences in trait expression between P-limited/co-limited plots and N-limited plots (Hypothesis 1). In detail, in plots with P limitation/co-limitation, individuals of Lythrum salicaria had a significantly shorter flowering period in individuals (Table 2), Plantago lanceolata started to flower later, and flowering of Lotus uliginosus and Anthoxanthum odoratum peaked later (Table 2). In these plots, Anthoxanthum odoratum had fewer inflorescences per individual and Plantago lanceolata had longer seed stalks and the height of the seeds above the soil surface was greater (Table 2). No other sexual reproduction traits of those species with significant difference were found. Besides, there was no significant difference in sexual reproduction traits of Holcus lanatus, Juncus acutiflorus, or Lysimachia vulgaris between P-limited/co-limited plots and N-limited plots (Table 2). To test our second hypothesis, we analyzed interspecific differences in trait expression between characteristic species from P limitation/co-limitation and characteristic species from N limitation. In species from sites with P limitation/co-limitation, the first flowering date, peak of flowering, and last flowering date were all significantly earlier than in the species from N-limited plots (Table 3).
When further exploring the intraspecific and interspecific differences in first weighted then averaged sexual reproduction traits along the plant and soil N:P, several strong correlations were found (especially along soil N:P) (Fig. 3). In detail, flowering phenology, i.e. first flowering date (FFD), peak of flowering (PF), last flowering date (LFD), flowering period in the population (FPP), and flowering period in individuals (FPI) were all significantly negatively correlated with increasing plant and soil N:P (Fig. 3a, b, c, d, e, f, g, h, i, j), indicating an earlier flowering time, but shorter flowering period under P limitation. Moreover, seed stalk length with seed panicle (SLP), seed panicle length (SPL), and percentage of flowered individuals (PFI) were all positively correlated with increasing soil N:P (Fig. 3l, n, p), while no effect of plant N:P was found on those three traits. Furthermore, negative effects of plant N:P were found on number of flowers (inflorescences) per individual (NF(I)I) (Fig. 3q). Note that the trends mentioned above are based on the species pooled without accounting for differences in species occurrence between the plots (see Appendix 1 Table 4). In general, we conclude from Fig. 3 that inter-as well as intraspecific variation in trait expression is large and that species behave very different in their trait expression along the nutrient gradients we measured.
Intraspecific and interspecific sexual reproduction performance of the selected grassland species along the gradients in soil pH and soil moisture We additionally explored the intraspecific and interspecific differences in trait expression along the soil pH and soil moisture gradients. We included soil pH and soil moisture since these environmental factors were found to differ significantly between the plots with P limitation/co-limitation and the N-limited plots (Table 1).
Strong correlations were found between the two environmental factors (especially soil pH) and first weighted then averaged sexual reproduction traits of all the species in each plot (Appendix 4, Fig. 5  (LFD), flowering period in the population (FPP), and flowering period in individuals (FPI) were all significantly negatively correlated with increasing soil moisture (Appendix 4, Fig. 5b, d, f, h, j), indicating an earlier flowering time, but shorter flowering period under high moisture condition. However, the traits mentioned above were all positively correlated with increasing soil pH (Appendix 4 Fig. 5a, c, e, g, i). Moreover, seed stalk length with seed panicle (SLP), seed panicle length (SPL), and percentage of flowered individuals (PFI) were all negatively correlated with soil pH (Appendix 4 Fig. 5k, m, o), while no effect of soil moisture was found on those three traits. Furthermore, negative effects of soil moisture were found on number of flowers (inflorescences) per individual (NF(I)I) (Appendix 4 Fig. 5r). Again, the trends mentioned above are based on the species pooled without accounting for differences in species occurrence between the plots (see Appendix 1 Table 4). In general, we also conclude from Appendix 4 Fig. 5 that inter-as well as intraspecific variation in trait expression is large and that species behave very different in their trait expression along the soil pH and soil moisture we measured.

Discussion
In the field, we tracked and measured a list of plant traits that can be regarded as a set of proxies for the investment in sexual reproduction and dispersal capacity in ten plots that covered all three types of nutrient limitation (N limitation, N and P co-limitation, and P limitation) and were distributed over three grassland sites in the Netherlands. We did this to test whether there were significant intraspecific and interspecific variations in the expression of sexual reproduction traits of grassland species along an N:P gradient. In addition to plant N:P and soil N:P, we also measured and analysed soil pH and soil moisture in the 10 plots as two additional environmental parameters.

Nutrient conditions in the plant and in the soil of the investigated grassland plots
In this field survey, in which we used a method of measuring nutrient concentrations in above-ground biomass to show the general nutrient conditions for plant growth that has been applied in various studies (e.g. Fujita et al. 2014;Olde Venterink et al. 2003;Roeling et al. 2018;Verhoeven et al. 1996), we did indeed find a gradient of N:P stoichiometry, from N limitation (Middenduin) to N and P co-limitation (Koeiemeent), and P limitation (Laegieskamp). However, because there were fewer P-limited plots than N-limited plots, and because of the proximity and similar plant growth conditions of Koeiemeent (N and P co-limitation) and Laegieskamp (P limitation), we mainly focused on comparing differences between N limitation and P limitation/co-limitation. In addition, soil nutrients were measured in order to provide extra information to further indicate nutrient conditions. It turned out that the trends in soil N:P and plant N:P were quite similar among the selected plots (R 2 = 0.444) (see Online Resource Fig. S2i). Moreover, measurements of soil nutrients supported the exploration of relationships between nutrient conditions and sexual reproduction traits, e.g. seed stalk and seed panicle length, and percentage of flowered individuals were all more correlated with soil N:P rather than plant N:P in this survey (Fig. 3l, n, p). Species invest less in sexual reproduction traits under P limitation/co-limitation than they do under N limitation (intraspecific differences) In the present survey, unfortunately, none of the species was present in all plots and very few were present in all three areas, which made it difficult to study the sexual reproduction performance of the selected grassland species intraspecifically, along the full N:P gradient from N limitation to P limitation/co-limitation. Nevertheless, our results generally support our first hypothesis that under P limitation/co-limitation, species invest less in sexual reproduction than they do under N limitation.
Seven out of the 11 species occurred in sufficient plots to enable an analysis of intraspecific variation (intraspecific differences) of sexual reproduction trait expression between N limitation and P limitation/colimitation, but significant differences were found for only four species: Lythrum salicaria, Plantago lanceolata, Lotus uliginosus, and Anthoxanthum odoratum. In general, compared with species occurring under N limitation, individuals occurring under P limitation/co-limitation showed lower investments in sexual reproduction (i.e., shorter flowering period in individuals (Lythrum salicaria), later flowering starting time (Plantago lanceolata), later peak of flowering (Lotus uliginosus and Anthoxanthum odoratum), and a lower number of inflorescences per individual (Anthoxanthum odoratum) ( Table 2). These species responses are in line with the finding of a greenhouse experiment with Holcus lanatus, a grass species that showed lower expression of several sexual reproduction traits at low P supply (Wang et al. 2019). Surprisingly, however, in the present field survey, Holcus lanatus showed no intraspecific differences in sexual reproduction performance (between the co-limited plot (KM1) and the N-limited plots). The possible reason might be the absence of this species in the most strongly P-limited area in our study (Laegieskamp). Other research has also shown that deficiency of P leads to a strong reduction in the number of flowers produced in a range of wetland species (Brouwer et al. 2001). The low investment in sexual reproduction under conditions of high N:P in our study further supports the idea that plant species can apply a strategy of economizing their P expenditure in response to P deficiency by investing less in sexual reproduction (Fujita et al. 2014;Wang et al. 2019), since classical allocation studies (Fenner 1986;van Andel and Vera 1977) have shown that sexual reproduction generally requires a higher percentage of P than of N and other elements plant acquired.
However, unexpectedly and contrary to our hypothesis, Plantago lanceolata in the plots with P limitation/co-limitation had higher and longer seed stalks (SHP and SLP) than it had in N-limited plots (Table 2), whereas in our greenhouse experiment we found that P limitation restricted the development of the seed stalk length and height of Holcus lanatus (Wang et al. 2019). To our knowledge, although seed release height is the most important plant-controlled dispersal parameter for grassland plants (Soons et al. 2004), studies of the influence of nutrient limitation on seed stalk length and height are lacking, possibly because such measurements are labor-intensive.
Moreover, in our study, under different nutrient conditions, the responses in sexual reproduction traits differed between species. With regard to flowering time, for instance, Lotus uliginosus and Anthoxanthum odoratum responded to P limitation/co-limitation in terms of peak of flowering time (the date with the most flowers or inflorescences), whereas the only significant flowering time response in Plantago lanceolata was a significantly postponed first flowering date (the date on which the first inflorescence occurred in that community) ( Table 2). The differences in species responses in our field study coupled with the lack of systematic experimental studies on the relationship between trait expression and nutrient supply and limitation in a wide range of species make it hard to draw general conclusions about the intraspecific plasticity of sexual reproduction traits along nutrient gradients in the field. Interspecific differences of sexual reproduction traits of the selected grassland species along an N:P gradient The results of our field survey also allow for the analysis of how trait expressions may change along an N:P gradient at interspecific level. With regard to flowering phenology, in characteristic species of environments with P limitation/ co-limitation, the first flowering dates (FFD), peaks of flowering (PF), and the last flowering dates (LFD) were all earlier than those of characteristic species of N-limited environments (Table 3; Fig. 3a, b, c, d, e, f). This is in contrast with previous research that reported that species persisting in high N:P had later flowering starting time (Fujita et al. 2014). Moreover, seed stalk length (SLP) and seed panicle length (SPL) were longer in P limitation/co-limitation (Fig. 3l, n), which were opposite to the result of Wang et al. (2019) that Holcus lanatus had shorter seed stalks and seed panicles under P limitation. A possible reason for the earlier flowering time and longer stalk length and seed panicle length we found under conditions of P limitation/co-limitation might be the influence of soil pH, which is known to have a significant effect on plant growth as well as on nutrient availability to plants (Bernal and McGrath 1994;Evans and Conway 1980). In our survey, N-and P-limited sites also differed significantly in soil pH, with low pH coinciding with significantly higher soil N and significantly lower soil P (Table 1, Appendix 2 Table 5). Furthermore, soil pH was significantly correlated with sexual reproduction performance of the grassland species: in acidic soils, flowering time was earlier, and seed stalks and panicles were longer (Appendix 4 Fig. 5). This opposite but possibly stronger influence of soil pH on investment in sexual reproduction than that of plant and soil N:P may have confounded the effect of high N:P on the investment in sexual reproduction reported by Fujita et al. (2014) and Wang et al. (2019), even though-as noted aboveat the intraspecific level, high N:P levels restricted flowering time significantly. Interestingly we found an earlier flowering time under higher soil moisture conditions, with P-limited/co-limited and lower pH conditions (Appendix 4 Fig. 5b, d, f). However, a high soil moisture usually coincides with lower soil temperatures (Miralles et al. 2012), and lower ambient temperatures usually lead to later flowering time (Jagadish et al. 2016). So, we probably reject the possibility that high soil moisture was a reason for earlier flowering times in these plots (Appendix Tables 4 and  5).
On the other hand, flowering period in the population (FPP) and in individuals (FPI), as well as number of flowers (inflorescences) per individual (NF(I)I) were restricted by P limitation/co-limitation (Fig. 3g, h, i, j, q). The restriction of P limitation/co-limitation on flowering period is consistent with the conclusion of Fujita et al. (2014) that species persisting in P-limited communities had shorter flowering period than those of N limitation/co-limitation. This phenomenon also agrees with the result of a study on Holcus lanatus, intraspecifically indicating that flowering period was shorter at both the population level and the individual level under P limitation (Wang et al. 2019). In addition, higher soil moisture was found to correlate with shorter flowering period (Appendix 4 Fig. 5h, j). Moreover, when soil moisture was higher, the number of flowers (inflorescences) per individual was lower (Appendix 4 Fig. 5r). The shorter flowering period, and less flowers (inflorescences) under higher soil moisture may indicate a restrictive influence of higher soil moisture on the investment in sexual reproduction of grassland species, given that soil moisture has been shown to influence various plant physiological processes via ambient temperature as also mentioned above (Hood 2001;Miralles et al. 2012).

Conclusions
In line with our first hypothesis, there was significant intraspecific variation in sexual reproduction traits along the N:P gradient we studied, with expression of reproduction traits being generally lower under P limitation. However, which specific trait showed plastic responses along the gradient differed per plant species, making it difficult to interpret these results as providing unequivocal support for the notion that plants adopt a general strategy of reducing their investment in sexual reproduction, thereby economizing on their use of P to be able to persist under P-limited conditions.
In line with our second hypothesis, there was significant interspecific variation in sexual reproduction traits along the N:P gradient. P limitation/co-limitation led to shorter flowering period and lower number of flowers (inflorescences) per individual. However, earlier flowering time, and longer seed stalk and seed panicle length were also found under P limitation/colimitation. The reason of the promotion of P limitation/co-limitation on earlier flowering time and longer seed stalk and seed panicle length could probably be the confounding effect of soil pH. Moreover, it must be noted that the interpretation of trait expression patterns between plots may have been hampered by the limited number (only 11) of species included in the comparison. We deliberately restricted the number of species studied because it is extremely arduous to monitor these traits in a large number of individuals in a field situation.
However, this research cannot disentangle the confounding effects of soil pH and soil moisture from soil N:P ratio, due to the coincidence of the P-limited/ co-limited environments and low soil pH and higher soil moisture. Therefore, we recommend that more indepth field surveys will be carried out in future on a wider range of species and longer nutrient gradients. Such studies are currently rare because of the practical challenges associated with large-scale field monitoring of phenological traits. Complementing field studies with controlled fertilization experiments may alleviate those challenges and could shed more light on the potential confounding effects of other environmental variables covarying with nutrient availability in the field. Another recommendation for future research comes from the observation that our 10 plots are located in only three different nature areas. It could be that sexual reproduction traits might have been evolved differently in these three areas over hundreds of years or longer, which might imply that differences of sexual reproduction traits could be attributed to the genetic differentiation, especially at the intraspecific level. Testing such a hypothesis requires common garden experiments in which the studied species from the three populations should be sowed or planted in an attempt to detect the possibility of confounding effects of genetic differentiation.