text:
Anatomical structure interpretation of the effect of soil environment on fine root function

Fine root anatomy plays an important role in understanding the relationship between fine root function and soil environment. However, in different soil environments, the variation of fine root anatomical structure in different root sequences is not well studied. We measured the soil conditions and anatomical structure characteristics (root diameter, cortical tissue, vascular tissue and xylem) of fine roots of Cupressus funebris in four experimental sites, and analyzed each level of fine roots separately. We link these data to understand the relationship between fine root anatomy and soil conditions. We found that the anatomical structure of fine roots is closely related to soil environmental factors. The fine roots of lower root order are mainly affected by soil nutrients. Among them, the cortical tissue of first-order fine roots was positively correlated with potassium and phosphorus, but negatively correlated with nitrogen, while second- and third-order fine roots was positively correlated with soil total potassium and negatively correlated with nitrogen and phosphorus. For the fine roots of high root order, the cortical tissue disappeared, and the secondary vascular tissue was mainly affected by soil moisture. In addition, we also found that the division of fine root functional groups is not fixed. On the one hand, the function of third-order fine roots will slip. For example, the decrease of soil moisture will promote the transformation of third-order fine roots into transport roots, and the reduction of nitrogen will promote the transformation of third-order fine roots into absorption roots to fix nitrogen. This transformation strategy can effectively prevent the restriction of soil nutrients on plant growth. On the other hand, with the change of habitat, the first- and second-order fine roots are still the absorbing root, and the fourth- and fifth-order fine roots are still the transport root, but the efficiency of absorption and transport will be affected. In conclusion, our findings emphasize the fine roots in different soil environment to show high levels of plasticity, shows that fine root anatomical structure changes may make plants, and reveals that the fine is just order of reaction and its mechanism in the soil environment.

Fine root anatomy plays an important role in understanding the relationship between fine root function and soil environment. However, in different soil environments, the variation of fine root anatomical structure in different root sequences is not well studied. We measured the soil conditions and anatomical structure characteristics (root diameter, cortical tissue, vascular tissue and xylem) of fine roots of Cupressus funebris in four experimental sites, and analyzed each level of fine roots separately. We link these data to understand the relationship between fine root anatomy and soil conditions. We found that the anatomical structure of fine roots is closely related to soil environmental factors. The fine roots of lower root order are mainly affected by soil nutrients. Among them, the cortical tissue of first-order fine roots was positively correlated with potassium and phosphorus, but negatively correlated with nitrogen, while second-and thirdorder fine roots was positively correlated with soil total potassium and negatively correlated with nitrogen and phosphorus. For the fine roots of high root order, the cortical tissue disappeared, and the secondary vascular tissue was mainly affected by soil moisture. In addition, we also found that the division of fine root functional groups is not fixed. On the one hand, the function of third-order fine roots will slip. For example, the decrease of soil moisture will promote the transformation of third-order fine roots into transport roots, and the reduction of nitrogen will promote the transformation of third-order fine roots into absorption roots to fix nitrogen. This transformation strategy can effectively prevent the restriction of soil nutrients on plant growth. On the other hand, with the change of habitat, the first-and second-order fine roots are still the absorbing root, and the fourth-and fifth-order fine roots are still the transport root, but the efficiency of absorption and transport will be affected. In conclusion, our findings emphasize the fine roots in different soil environment to show high levels of plasticity, shows that fine root anatomical structure changes may make plants, and reveals that the fine is just order of reaction and its mechanism in the soil environment. KEYWORDS fine root, anatomical structure, root function, root order, soil environment, Cupressus funebris Introduction Fine roots are the most physiologically active plant components of the root system (Cuneo et al., 2020). With short life span and fast turnover rate (Hu et al., 2021), fine roots are mainly responsible for absorbing a large amount of water and nutrients and carrying out carbon (C) and nitrogen (N) cycling (Zeng et al., 2020;Endo et al., 2021;Huaraca et al., 2021), which ensuring an adequate supply of water and nutrients for plant photosynthesis, growth, and maintenance (Freschet and Roumet, 2017). Fine root anatomical structure is the direct expression of root growth and development, and there is a strong relationship between anatomical structure and physiological function Zhou et al., 2022). The anatomical structure of fine roots is mainly composed of cortex and stele (Guo et al., 2008), which reflects two main functions of fine roots: absorption and transport Mao et al., 2018;Wang et al., 2019). In recent years, root functional ecology and physiological ecology have attracted much attention (Freschet et al., 2021). However, despite the recognition of functional differences, few studies have involved anatomical structures, especially the anatomical structures of fine roots in different soil conditions. Due to the difficulty of sampling and the complexity of research, the research on fine root anatomical structure and its relationship with environment and related physiological activities is not sufficient.
Root order is vital for understanding the root function in fine root system (Guo et al., 2008). Traditionally, fine roots generally consist of 2-5 root orders (Freschet and Roumet, 2017). Specifically, roots farthest from the main root axis of the root system, being root tips with no more branches, were defined as the 1-order fine root; the parent root of the 1-order fine root was defined as the 2-order fine root; and so on, until the 5-order fine root was reached. In addition, roots borne on a higher order fine root without branches were also classified as a 1-order fine root (He et al., 2022). Many studies have emphasized the differences of anatomical structures of fine roots in different branch orders (Liu et al., 2019). From first-to fifth-order fine roots, the cortex gradually decreases (or disappears), while the stele gradually increases (Guo et al., 2008;Endo et al., 2021). Later, the researchers found that the functional characteristics of fine roots with different branch order can be directly expressed by the internal anatomical structure of fine roots (Guo et al., 2008). Specifically, transport fine roots are generally high order, have thicker diameter and greater ratio of stele to root diameter. In contrast, absorptive fine roots are low order, have thinner diameter and lower ratio of stele to root diameter . That is to say, from 1-to 5-order fine roots, the absorption capacity of fine roots gradually decreases, while the transport capacity gradually increases (McCormack et al., 2015). In fact, this change in the physiological function of advanced fine roots is related to secondary development , such as the development of cork periderm, endodermis (Wells and Eissenstat, 2002;Hishi, 2007). This is because with the increase of root order, fine roots require greater axial transport capacity. These changes in anatomical structure reduced the ability of roots to absorb water and nutrients, and increased the transport efficiency of higher roots .
Ecological factors can affect the anatomical structure of fine roots Huaraca et al., 2021). Temperature can affect not only plant roots directly, but also fine roots indirectly by affecting soil nutrient conditions (Sorensen et al., 2008). For example, at low temperature, plants get more water and nutrients by increasing the diameter of fine roots and producing more absorbing roots (Zadworny et al., 2016), so as to improve the adaptability of fine root to environmental change. It is also worth noting that fine root growth also increases significantly at high temperature (Malhotra et al., 2020). Present day science indicates that the stimulation of photosynthesis at high temperature and the prolongation of growing season were beneficial to root growth , and the response mechanism of fine roots to climate change was further clarified. In addition, soil moisture also has a certain effect on the anatomical structure of fine roots. The absorption of water by the root first involves the radial water transport from the soil solution to the root column, and then flows upward into the shoot through the xylem catheter . The response of fine roots to drought involves the complex relationship between anatomical structure and function (Barrios-Masias et al., 2015;Reingwirtz et al., 2021). Studies have shown that plants in arid environments will carry out drought avoidance strategies . Under drought conditions, the roots need to absorb enough soil water to meet the transpiration needs of the canopy (Cuneo et al., 2020). At this time, the fine root enhances the absorptive capacity by reducing the diameter (Nikolova et al., 2020). The results showed that the root system adapted to soil water stress changed the anatomical structure, so as to achieve the best growth state.
Further, the fine root anatomical structure also changes with the local soil environment (Hogan et al., 2020). Soil bulk density is also likely to be an important driver of fine root trait variation , it usually has a positive effect on the fine root diameter (Clark et al., 2003). Plants with higher fine root diameter have greater root elongation performance in denser soils , whereas the contrary may be true in less dense soils . This is most likely because a smaller soil bulk density provides a lower mechanical impedance, which promotes root dispersal (Cao et al., 2021).
At the same time, it is found that larger soil pores and fine root density can create more preferential flow paths for stem flow and natural precipitation (Johnson and Lehmann, 2006), promoting plants to absorb more water and nutrients. Additionally, about soil nutrients (nitrogen, phosphorus, potassium) and anatomical structure of fine roots, lack of nutrients in the soil usually leads to root elongation and structural development (Hogan et al., 2020). Nitrogen (N) is arguably the most studied element in plant function and growth (Hogan et al., 2021). In the soil with high concentration of N, the metabolic activity of fine roots is stronger Frontiers in Plant Science 03 frontiersin.org (Bardgett et al., 2014), such as root elongation , nutrient uptake (Freschet and Roumet, 2017), proton release and active C exudation (Jones et al., 2009). With the increase of N content, the cortical thickness of the first three grades of fine roots decreased significantly (Chen et al., 2010), indicating that the increase of N content could promote the transport function of fine roots. In general, N deficient soil leads to increased lateral root elongation, while phosphorus (P) deficient soil leads to more branching roots (Giehl et al., 2013). On the one hand, the decrease of soil available phosphorus will reduce the diameter of fine roots, on the other hand, the surface area of fine roots will increase (Ushio et al., 2015). Therefore, fine roots can still absorb water and nutrients to the maximum extent. Several other micronutrients, such as calcium (Ca), magnesium (Mg), sulfur (S), and other trace metals, are required in small amounts by plants but rarely limit their physiological functions. For example, calcium, if they are in compounds that are phloeminsoluble or too large for phloem transport, since the amount of calcium in the soil has relatively little effect on plants (Hogan et al., 2021). In general, the root diameter decreased in areas with low soil nutrients (Hogan et al., 2020). These studies show that with the changes of soil environment, the distribution and morphological structure of the fine root system will change obviously (Hertel and Schöling, 2011;Moser et al., 2011;Girardin et al., 2013;Mao et al., 2015).
Taken together, these studies support the notion that fine roots can adjust a variety of structural and functional traits to adapt to different soil conditions, which is closely related to the changes of fine root anatomical structure and is also vital for the growth of trees and the sustainable development of forests. Considering that the variation of fine root anatomical structure will lead to the change of fine root function, we suppose the fine root function group will change in different soil environment. Meanwhile, although some scholars have studied the anatomical structure based on artificially changing some soil physical and chemical indexes, how do fine roots adapt to the long-term natural growth of soil by changing the anatomical structure (Nikolova et al., 2020)? It's not clear yet. In this case, it is crucial to determine how fine root anatomical structures of different root orders adapts to and responds to changes in soil environment and the relationship between fine root anatomical structure and function.
Our study is based on Cupressus funebris of several distribution areas including Suining, Mianyang, Deyang and Guangan, which represented the hilly areas in the middle of Sichuan. We excavated roots by the soil monolith method . In this study, we attempted to find out the changing trend and response mechanism of anatomical structure under different soil environmental factors by measuring the anatomical structure of first-to fifth-order fine roots between different sites and the soil physical and chemical properties of the corresponding sites. Based on this, we propose two hypotheses: (1) the fine root anatomical structure of C. funebris in four sites showed plasticity, the response of fine root anatomical structure of different root orders to soil environmental factors is not consistent, (2) due to changes in anatomical structure, the division of fine root functional groups is not fixed.

Study site' description
We selected four sampling sites with large differences in soil environment in rural mountainous areas, namely, Jingyang District of Deyang city (DY), Anju District of Suining City (SN), Yanting County of Mianyang City (MY) and Huaying City of Guang'an City (GA) to represent shallow hill area, middle hill area, deep (high) hill area and low mountain area successively ( Figure 1). These four sites are all ecological groups of singlespecies C. funebris plantations, located in the subtropical humid monsoon climate zone. The forest topsoil thickness is about 30 cm, and the age of C. funebris is 25-30 years old. We choose this age of C. funebris because at this time the growth of C. funebris is relatively stable, which is more conducive to the determination of anatomical structure indicators. The dominant species of understory plants were Myrsine africana, Vitex negundo and Coriaria nepalensis. The detailed nature of the test site is shown in Table 1.

Experimental design, root sampling and anatomical structure analysis
Our experimental study encompasses four sites. According to the preliminary investigation (stand density, tree height, diameter at breast height (DBH), natural conditions, etc.), in each site, we randomly selected three standard plots with light disturbance, well-developed natural vegetation, canopy density of about 0.8 and similar vertical space as standard plots (He et al., 2022).Then, 3 standard plots × 4 sites, for a total of 12 established in this study, with a standard plot area of 600 m 2 (20 m × 30 m).
The fine roots sampling took place in July 2017. According to DBH and tree height, we selected 3 relatively scattered trees in each plot as standard trees. Then, we selected four directions in each standard wood (there was no obvious direction in the fine root of C. funebris wood after the previous study), and determined that it was the main root of the plant according to the direction of the root, dug carefully along the main root, and collected a complete root segment containing 1-to 5-order fine roots born on the main root. In each sample sites, we selected 18 1-order fine roots, 18 2-order fine roots, 15 3-order fine roots, 6 4-order fine roots, and 6 5-order fine roots, and four sites have a total of 756 roots. After collection, the fine roots samples were stored in a refrigerator at low temperature (0-4°C) and transported to the laboratory as soon as possible.
In the laboratory, we carefully separate the fine roots from the soil and gently removed soil particles and debris that cling to the roots. The root samples were washed with deionized water at room Frontiers in Plant Science 04 frontiersin.org temperature and the fine roots were stored as anatomical samples in formalin-aceto-alcohol (FAA) fixation solution. If the root segment was too long, we processed it to about 1 cm to facilitate subsequent experimental operations. Then we use the paraffin technique to make fine root cross sections, and all root cross sections were stained with safranin-fast green. After gradient alcohol dehydration (alcohol concentrations of 70, 85, 95 and 100%), safranin-fast green staining, wax-impregnation (using xylene as a drying solvent at 45°C in an oven for 20 h), and paraffin embedding, fine roots were cut into thin sections, and dewaxed to achieve transparency. Finally, 8 μm thick sections were prepared for analysis (Liu et al., 2019). For each root cross section, eight sampling points on each slice with complete root morphology were selected and photographed under a polarizing microscope (BX51; Olympus Corporation, Tokyo, Japan). The root diameter (RD) (mm), vascular bundle diameter (VBD) (mm), cortex thickness (CT) (mm), xylem area (XA; cm 2 ), number of vessels (NV) were measured from one root cross section order sample. In addition, the stele: root diameter ratio (VBD/RD) is calculated to represent the ratio of root diameter used for resource transportation. And we calculated cortex ratio (C/R; cortex area / root cross section area) from the measured (Natsuko et al., 2019).

Soil sampling, physical, and chemical properties analysis
First of all, the button thermometer (DS1921G) was covered with a sealed bag and buried in the soil to continuously measure the soil temperature, set to record every 2 h. Then, we used the ring-knife method (inner diameter 5 cm) to collect soil samples 50 cm from the trunk of each standard tree for the determination of soil bulk density (SBD) and soil porosity (SP). Finally, after removing the litter and humus layer on the soil surface,  we collected three topsoil samples with a 5 cm diameter drill at a depth of 0-20 cm at the fine root sampling point for the determination of soil chemical properties. This depth was chosen since fine roots are mainly confined to the upper soil, and the relationship between fine roots and physical and chemical properties of soil should be stronger than that of deep soil (He et al., 2022). After collection, all soil samples were stored in a refrigerator at low temperature (0-4°C) and transported to the laboratory as soon as possible for further analysis.
In the laboratory, we mixed all soil samples were mixed into one composite sample per tree to determine the physical and chemical properties of soil (He et al., 2022). First, the physical properties of the soil: SBD and SP were determined by the ringknife method (Yokobe et al., 2021); SW was determined by determining water loss (Yokobe et al., 2021). Then, remove impurities (such as rocks, plant litter) and determine the chemical properties of the soil: soil total nitrogen (STN) was determined by the semi-micro Kjeldahl digestion method (Sinha et al., 1986); soil total phosphorus (STP) was determined by the sodium-hydroxide fusion molybdenum antimony anti-coloring method (Fisher et al., 1998;Hogan et al., 2021); soil total potassium (STK) was determined by the sodium-hydroxide melting-flame atomic absorption spectrophotometry method (Hogan et al., 2021); soil alkali-hydrolyzed nitrogen (SAN) was determined by the alkalisolution diffusion method (Hogan et al., 2021); soil available phosphorus (SAP) was determined by the sodium bicarbonate extraction-molybdenum-antimony colorimetric method (Cao et al., 2021); soil available potassium (SAK) was determined by the neutral ammonium acetate extraction-atomic absorption spectrophotometry determination method (Hogan et al., 2021); soil organic carbon (SOC) was determined by the potassium dichromate acid oxidation-external heating method (Curtin and Smillie, 1986).

Data analysis
First, the Shapiro-Wilk test and Levene's test were used to determine the normality and homogeneity of variance for each group of variables, respectively. And use Excel software (Version 2019, Microsoft Office, United States) to manage, sort and calculate the data (He et al., 2022).
Second, the SPSS software (Version 25.0, SPSS Inc., Chicago, IL, United States) was used to analyze the anatomical characteristics of fine roots of different root orders (RD, VBD, VBD/RD, CT, C/R, XA, NV) and soil physical and chemical properties (SOC, STN, STP, STK, SAN, SAP, SAK, SW, SBD and SP). The differences in fine roots anatomical structure, soil physical and chemical properties, and root orders among the four test sites were analyzed by one-way ANOVA, followed by Duncan's multiple comparison method (p < 0.05). Two-way ANOVA was used to analyze the effects of soil factors, root order and their interactions on fine root anatomical structure and nutrients.
Finally, the redundancy analyses (RDA) (Li et al., 2020) was performed to analyze the relationships in fine roots anatomical structure, soil environmental factors, and nutrients content by using Canoco software (version 5.0).

Results
Characteristics of fine root anatomical structure in the different city sites First, we analyze whether different soil conditions have an impact on the anatomical structure of each level of fine roots. Our results show that both the four experimental sites and root order have a significant influence on the anatomical structure of fine roots (p < 0.01, Table 2). However, the interaction between them has different influences on them, specifically, it has extremely significant influence on VBD and VBD/RD (p < 0.01, Table 3), significant influence on XA and NV (p < 0.05, Table 2), but no influence on RD, CT and C/R (p > 0.05, Table 2). That is to say, the interaction between sites and root order mainly affects the stele, but not the cortex. So the interaction between them mainly affects the transport capacity of fine roots, but has little effect on the absorption function.
In the four experimental sites, RD, VBD, VBD/RD, CT, C/R, XY and VE gradually increased with the increase of root order. For the first-order fine roots of C. funebris, there are significant differences in RD, VBD, VBD/RD and C/R between MY, DY and SN, GA (p < 0.05, Table 3). For the second-order fine roots, RD, VBD/RD, C/R and NV are significantly different between MY, DY and SN, GA (p < 0.05, Table 3). For the third-order fine roots, RD, VBD/RD and C/R are significantly different between MY, DY and SN, GA (p < 0.05, Table 3). For the fourth-order fine roots, VBD and VBD/RD were significantly different between MY, DY and SN and GA, and the contents of XA and NV in SN were significantly higher than those in the other three regions (p < 0.05, Table 3). For the fifth-order fine roots, RD is significantly different between MY and DY, SN, GA, VBD is significantly different between MY and GA, XA and NV are significantly different between MY, DY and SN, GA (p < 0.05, Table 3). The differences in anatomical structure of fine root sequences between the four sites indicated that fine roots had heterogeneous plasticity.

Correlation between fine root anatomical structure and soil environmental factors
Secondly, the physical and chemical properties of soil in the four experimental sites are different after our analysis. PCA analysis showed that the soil nutrients were the highest in DY, followed by MY and SN, and the lowest in GA (He et al., 2022). The details is shown in Supplementary Table S1. For the first-order fine roots, vascular bundle and xylem are positively correlated with N and SOC, negatively correlated with K and SW, while the cortex is just the opposite. And the diameter was positively correlated with K, SW and SBD, and negatively correlated with N, P and organic carbon (Figure 2A). For   TABLE 2 Variance analysis of the influence of site and roots order on the anatomical structure of fine roots.  second-order fine roots, vascular bundles and xylem were positively correlated with N and SOC, and negatively correlated with K, SBD and SW, while the cortical tissue was on the contrary. And the diameter is negatively correlated with SW, SBD and SP ( Figure 2B). For third-order fine roots, vascular bundles, the tissue and xylem were positively correlated with N, P, SOC and SBD, and negatively correlated with SW. The cortical tissue is just the opposite. The root diameter was positively correlated with SW and negatively correlated with N and SOC ( Figure 2C). The fifth-order fine roots are very similar to the fourth-order fine roots. There was no cortical tissue in all of them. Vascular bundle tissue and xylem were positively correlated with N, P, K, SOC and SBD, and negatively correlated with SW. The difference is that the diameter of fourth-order fine roots are also positively correlated with N, P, and the diameter of fifth-order fine roots are positively correlated with K and SW ( Figures 2D,E). These results show that the anatomical structure of fine roots at all orders is closely related to soil factors, showing specific adaptation strategies of fine roots. Finally, the importance of soil environmental factors to the anatomical structure of 1-5 order fine roots was ranked. The anatomical structure of first-order fine roots of C. funebris was mainly affected by SAN, STK (p < 0.05), accounting for 53 and 8.6% of the total information of soil environmental factors, respectively. Second-order fine roots were most affected by SW, SAK and STP (p < 0.05), accounting for 42.6, 19.7 and 11.9%, respectively. Third-order fine roots were mainly affected by SW and STP (p < 0.05), accounting for 60.3 and 10.6%, respectively. SW had the greatest effect on fourth-to fifth-order fine roots (p < 0.05), accounting for 60 and 55.9%. From this, we can conclude that the fine roots of high-order roots are mainly affected by soil moisture, while the fine roots of low-order roots are greatly affected by soil nutrients (N, P, K).

Anatomical structure of fine roots under different soil conditions
In response to the change of soil environment, the anatomical structure of fine root shows plasticity in the trade-off between the acquisition and protection of soil resources (soil water and nutrients; Liu et al., 2019;Yang et al., 2019), which supports our first hypothesis. Specifically, the anatomical characteristics of roots vary according to root sequence , soil nutrients (Chen et al., 2010), soil bulk density  and soil biological community (Hishi, 2007).
Our results show that the anatomical structure of fine roots of C. funebris shows obvious geographical differences with the change of soil environment. Cortical thickness is one of the important indexes to measure the ability of root absorption and radial transport (Guo et al., 2008). In this study, the areas with high STN and SAN (such as DY and MY, Supplementary Table S1) exerted thicker cortex, but smaller vascular bundle diameter and xylem area (Table 3). This is due to the decrease of cortical cell diameter or the number of cortical cell layers under high N condition (Oscar et al., 2001). It is also worth noting that, contrary to the response of cortical tissue, the diameter or cross-sectional area of fine roots increased under the condition of high N (Stefan et al., 2001). The fine root vascular bundle diameter, xylem tissue and vessel number reflect the root transport capacity. The higher the value, the stronger the fine root transport capacity and the stronger the ability to resist adverse environmental factors (Guo et al., 2008;Jia et al., 2013). Therefore, we found that in the environment with relatively high soil nitrogen content, the transport capacity of fine root was stronger than that in the environment with relatively low nitrogen content. At the same time, there is a correlation between the anatomical structure and physiological variation of fine roots (Weemstra et al., 2016). Several studies have shown that fine roots can respond to stressful environments by altering their anatomical structure, including the development of larger vascular bundles in saturated soils (Yamauchi et al., 2021) and the formation of aerenchyma (Kotula et al., 2015). This means that the development of cortex tissue of C. funebris may be relatively better in the environment of soil nutrient deficiency, while the development of vascular tissue of C. funebris may be relatively better in the environment of soil nutrient abundance.
Diameter is the most important morphological feature of root system with a secondary predictor of root anatomical structure variation (Gu et al., 2014). According to the growth law of plants, the thickness of diameter is determined by both the thickness of cortex and the diameter of vascular bundles. Most studies have shown that there is a significant positive correlation between fine root anatomical structure (cortex thickness, vascular bundle diameter) and root diameter in tropical, subtropical and temperate forests Liu et al., 2018). In this study, fine roots exerted larger diameter and thicker cortex in SN and GA, especially in 1-to 2-order fine roots (Table 3). We speculate that the reason may be related to the higher soil moisture content in SN and GA (Supplementary Table S1). Due to the root anatomical and physiological activities are closely related (McCormack et al., 2015), fine root morphological changes could be the result of water demand and physical changes (Kong et al., 2017). Studies have shown that root diameter decreases organization to reduce the main reason is that the root cortex (Kong et al., 2017(Kong et al., , 2019. Therefore, in areas with higher soil moisture content (such as SN and GA), the metabolic activity of fine roots is more vigorous, and to ensure the rapid acquisition and transportation of water and nutrients, the thickness of the cortex is also larger, and the fine roots are thicker. In addition, a thicker diameter may be built to maintain its absorptive capacity for a longer time (Zadworny et al., 2016), and the long-term investment in organizational construction and maintenance is relatively low (Valverdea-Barrantes et al., 2020).
At the same time, soil moisture was the main soil environmental factor affecting the anatomical structure of secondto fifth-order fine roots (Figure 2). Vascular tissue was negatively Frontiers in Plant Science 08 frontiersin.org correlated with soil moisture, while cortical tissue was positively correlated with soil moisture (Figure 2). The function of the xylem conduit is to collect and transport water and minerals from the root to the above ground . On the one hand, the water in the soil enters the roots through xylem "pipes" (Warren et al., 2015). And through the availability of soil water, root growth and microbial activity can be regulated, which indirectly affects the C cycle of underground soil (Curiel et al., 2003;Weltzin et al., 2009). On the other hand, drought-tolerant plants promote the formation of periderm in fine roots in higher root order (Blokhina et al., 2003). C. funebris has strong drought tolerance. In the environment with relatively little soil moisture (DY and MY), fine roots can promote the thickening of periderm thickness to reduce water loss and promote the development of xylem structure, which can transport water needed by various organs more effectively. Thus, the growth and development of root system is affected by environmental conditions have plasticity, thus to better face the changes in the environment.

Response of fine root order function in different soil environments
Root order rather than diameter has been shown to be a better representative of root function (Weemstra et al., 2016).
The differences in root anatomy among different root sequences enable them to perform different ecological functions, such as nutrient acquisition and transport (Yu et al., 2022). Usually, the five-level root order can be divided into two groups: the root whose main function is to absorb resources and the root that serves functions other than absorption (Guo et al., 2008). Absorptive fine roots are mainly involved in acquiring and absorbing soil resources, while transportable fine roots mainly play the functions of structure and transportation, and have some additional storage capacity (McCormack et al., 2015). The development of root function from absorption to transport may be at the cost of reducing the thickness of cortical tissue, indicating the difference of absorption function from high to low and transport function from low to high along the root order. Changes in fine root function usually represent different strategies or adaptations of plants to various environmental conditions (de la Riva et al., 2021).
Although the root order is different, plants have been shown to respond to the soil environment. With regard to the classification of fine roots, generally speaking, the first-to thirdorder fine roots of plants are low-order roots (absorbing roots) and fourth-to fifth-order fine roots are high-grade roots (transport roots; Wang et al., 2019). It has been realized that this classification method is too general to take into account the . RD, means fine root diameter, VBD, means vascular bundle diameter, VBD/RD, means stele: root diameter ratio, CT, means cortex thickness, C/R, means cortex ratio, XA, means xylem area, NV, means number of vessels; SOC, means organic carbon, STN, means Total nitrogen, STP, stands for total phosphorus, STK, stands for total potassium, SAN, stands for alkali hydrolyzed nitrogen, SAP, stands for available phosphorus, SAK, stands for available potassium; SW, stands for soil moisture, SBD, stands for soil bulk density, and SP, stands for soil porosity. The letters in this figure indicated the same meaning as the chart below.
Frontiers in Plant Science 09 frontiersin.org impact of soil environmental changes on the anatomical structure of fine roots, and cannot explain the changes in fine root functional groups (Freschet and Roumet, 2017). Recent studies have further emphasized that the structural and morphological characteristics of different root sequences can change independently of each other, making woody plants adapt to different environments .
With the change of soil environment, the function of different root orders is not fixed. For example, in the site with the highest STN (that is DY), the fine root diameter and the vascular bundle diameter of 3-order fine roots was high (Table 3). On the one hand, N was previously thought to play a similar core role in the characteristic tissues of roots, supporting metabolic activities (including nutrient and water transport, enzyme function and mycorrhizal; Ma et al., 2018); On the other hand, soil moisture is vital for the diffusion of N to the root surface (Warren et al., 2015). Therefore, we speculate that due to the low soil moisture content in DY (Supplementary Table S1), fine roots may not be able to transport enough N to the aboveground parts of plants through diffusion. Since, by increasing the diameter of fine roots and vascular bundles,3-order fine roots can also obtain enough N. And there are studies that show that as soil N increases, fine roots transport function increases by decreasing the thickness of cortex (Chen et al., 2010). Therefore, in order to fix nitrogen and absorb more water in arid soil, we analyzed that the 3-order fine roots have a tendency to evolve into transport roots. However, we observed that in areas with low N content (such as SN), the thickness of the 3-order fine roots cortex were significantly higher than that of the other three experimental areas (Table 3). In other words, in the soil with low N content, the 3-order fine roots of the plant may perform the function of absorption. The transformation of the function of 3-order fine roots under different soil conditions may be driven by environmental pressure to evolve fine roots that can adapt to different soils, which is a more competitive strategy for fine roots to compete for soil nutrients. Another key factor behind the occurrence of functional slippage in 3-order fine roots may be that the stele and cortical tissue are allometric (Valverdea-Barrantes, 2022). It was found that the ratio of the cortex (occupied root cross-sectional area) to the stele (occupied root cross-sectional area) increased as the root diameter increased will achieve a better balance between nutrient absorption and transport (Zhou et al., 2021). The allometric relationship followed by the stellars and cortex helps us to better understand the strategies for the acquisition and conservation of thirdorder fine root resources and their adaptation to different environments.
In addition, we analyzed that in our four different sites, 1-and 2-order fine roots showed intact cortex and low stele proportion, so we defined them as absorbent roots. Similarly, 4-and 5-order fine roots showed no cortex and obvious secondary development, so we defined them as transport roots (Guo et al., 2008). Of course, we also wonder whether the soil environment changes greatly, will the function of the 2-order and 4-order fine roots also change? For example, from tropical rain forest to desert or north forest. The root economic spectrum (RES) advocates a trade-off between maximizing access to resources and improving productivity, or maximizing protection of resources and prolonging life (Carmona et al., 2021;de la Riva et al., 2021). Studies have shown that the root life with high maximum absorption rate is shorter (Weemstra et al., 2016). From the perspective of cost-effectiveness, it is difficult to explain the mismatch between the fact that all 1-to 4-order fine roots are absorbing roots or 2-to 5-order fine roots are transporting roots and the conservation of root phylogeny. Therefore, 2-order fine roots will not become transport roots, and 4-order fine roots will not become absorbent roots. If the function of the 2-order and 4-order fine roots is also changed, then the absorption and transport functions of the plant will lose coordination. However, differences in the anatomical structure of fine roots at different levels affect the efficiency of absorption of nutrients and water in fine roots (1-and 2-order) and transportation of nutrients and water in fine roots (4-and 5-order).
These results further support the view that the responses of fine root anatomical structures with different root orders to soil environmental factors may not be consistent. This structure of the root system is helpful for the lower roots to more easily acquire and effectively utilize the resources in the soil environment, and the higher roots to more easily and efficiently transport water and nutrients to all parts of the plant, so as to maintain the demand for water and nutrients and normal growth of the plant.

Conclusion
We used several typical indexes (diameter, cortical tissue, vascular tissue and xylem) to evaluate the variation of fine root anatomical structure and the division of fine root functional groups with different root order under different soil conditions. Our conclusion is that soil environmental factors may influence the anatomical structure of fine roots with different root orders. First of all, for higher order roots, soil moisture (SW) was the most important factor; for lower order roots, soil nutrients (N, P, K) were more important. Then, our study fine roots are plastic in their response to local environmental conditions during development: low soil moisture and nutrient content promoted the development of fine root cortex, while high soil moisture and nutrient content promoted the development of vascular bundles. Finally, there was obvious heterogeneity in anatomical structure and function of fine roots with different root sequences. As fine roots adapt to different soil environments, the division of fine root functional groups may not fixed. Although different soil environments did not change the functional group characteristics of 1-2 order fine roots as absorption roots and 4-5 order fine Frontiers in Plant Science 10 frontiersin.org roots as transport roots, their absorption and transport efficiency would be affected. The exception is the 3-order root, which we define as a transitional root because it may be biased in favor of absorbing or transporting roots under different circumstances. In short, fine roots respond to changes in soil environment by changing different anatomical structures, and their corresponding functions will also change, which is the adaptive strategy of fine roots. The next step would continue to study the effects of soil environment on fine roots. For example, the nutrients and water physiology of fine roots are worth studying.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.


s2fieldsofstudy:
Agricultural And Food Sciences
Environmental Science

provenance: AgriculturalAndFoodSciences,2022-2022-Agricultural And Food Sciences-2022.jsonl.gz.gz

text:
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.