*Article* **Climatic Factors Determine the Distribution Patterns of Leaf Nutrient Traits at Large Scales**

**Xianxian Wang <sup>1</sup> , Jiangfeng Wang <sup>1</sup> , Liuyang Zhang <sup>2</sup> , Chengyu Lv <sup>3</sup> , Longlong Liu <sup>4</sup> , Huixin Zhao 1,\* and Jie Gao 1,5,\***


**Abstract:** Leaf nutrient content and its stoichiometric relationships (N/P ratio) are essential for photosynthesis and plant growth and development. Previous studies on leaf nutrient-related functional traits have mainly focused on the species level and regional scale, but fewer studies have investigated the distribution patterns of the leaf N and P contents (LN, LP) and N/P ratios (N/P) in communities and their controlling factors at a large scale; therefore, we used LN, LP, and N/P data at 69 sites from 818 forests in China. The results showed significant differences (*p* < 0.05) in the LN, LP, and N/P at different life forms (tree, shrub, and herb). Neither LN, LP, nor N/P ratios showed significant patterns of latitudinal variation. With the increase in temperature and rainfall, the LN, LP, and leaf nutrient contents increased significantly (*p* < 0.001). Across life forms, LN at different life forms varied significantly and was positively correlated with soil P content (*p* < 0.001). The explanatory degree of climatic factors in shaping the spatial variation patterns of LN and N/P was higher than that of the soil nutrient factors, and the spatial variation patterns of the leaf nutrient traits of different life forms were shaped by the synergistic effects of climatic factors and soil nutrient factors.

**Keywords:** leaf nutrients; functional traits; life forms; climate change; soil nutrients

### **1. Introduction**

The functional traits of plants reflect their survival strategies in response to environmental changes [1,2]. Nitrogen (N) and phosphorus (P) are important components of the basic structure of plant cells and their levels. The stoichiometric relationships between N and P drive photosynthesis, plant growth, reproduction, and other life processes [3]. Leaf nitrogen (LN), phosphorus (LP), and leaf N/P ratios (N/P) are key plant traits that influence the productivity of forest communities and regulate carbon cycling [4]. Numerous studies have found that climatic factors, soil nutrient factors, and community genealogical structure influence the plant functional traits by affecting plant metabolism [5]. Exploring the distribution patterns of functional traits related to leaf nutrients on a large scale has important ecological significance for quantifying the impact of environmental changes on plant functional traits [6]. However, the key drivers that determine the distribution patterns of these key functional traits remain elusive.

Climate factors influence the functional traits of plant leaves by controlling plant metabolism, growth, and development processes [7]. A global data-based study found that climate change significantly affected the leaf nutrient content [8]. With increasing temperature, the LN decreases significantly and LP increases significantly [9]. Higher LN

**Citation:** Wang, X.; Wang, J.; Zhang, L.; Lv, C.; Liu, L.; Zhao, H.; Gao, J. Climatic Factors Determine the Distribution Patterns of Leaf Nutrient Traits at Large Scales. *Plants* **2022**, *11*, 2171. https://doi.org/10.3390/ plants11162171

Academic Editors: Martina Pollastrini and Oleg Chertov

Received: 8 August 2022 Accepted: 18 August 2022 Published: 21 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

enhances photosynthesis under low-temperature stress, and when temperatures are higher, plants coordinate metabolic processes by reducing the LN and increasing LP [10]. Other studies have found that the N/P of plant leaves increased significantly with increasing temperature [11]. Numerous previous studies have shown that temperature is a key factor in determining changes in leaf nutrient traits, especially the mean annual temperature (MAT) [4,12,13].

Plant leaf nutrients are also influenced by precipitation factors [14]. Precipitation factors can alter soil nutrient availability, and when rainfall is low and soils are subject to drought stress, plants respond to water stress by changing functional traits [15]. Temperature and precipitation factors jointly determine the nutrient changes in plant leaves [16]. In specific forest communities, the key climatic factors that determine the leaf nutrients in different species remain elusive due to differences in biological adaptations and plant nutrient uptake strategies. Therefore, it is particularly important to quantify the effects of climatic factors in the leaves of different plant forms at a macroscopic scale.

Not only do climatic factors have an important impact on the nutrient profiles of plant leaves, but soil nutrient factors also play a key role in shaping the spatiotemporal pattern of leaf nutrient traits [17]. Climate affects the plant nutrient profiles by influencing soil nutrient redistribution [18]. Plant vessels absorb inorganic salts and water from the soil and transport them to the leaves, so the soil is the main source of nitrogen and phosphorus for plant leaves [19]. The nutrient profiles of leaves are limited by the nutrient availability of the soil [20]. Some studies have found that the N/P of Chinese plant leaves is significantly lower than the global average due to the low effective phosphorus content of Chinese soils [21]. It has been shown that the N and P content of plant leaves increases with increasing soil N content and soil pH [22], and that soil nutrient factors directly affect the nutrient profiles of leaves.

Forest communities are divided into tree, shrub, and herb levels according to species composition, structure, and production, with each plant having its life form, and vertical differentiation provides a good indication of the community's adaptation to environmental conditions [23]. The nutrient profiles of leaves vary significantly among different life forms due to differences in their survival strategies [24]. Comparisons of the leaf nutrient traits of different plant forms are mostly based on the species scale or local scale, and few studies have been conducted on a larger spatial scale based on the community scale. Community trait variation is more effective than species trait variation in predicting the plant response to environmental change [25]. However, due to interspecific competition and intraspecific struggle in a given forest community, there is a great theoretical risk in predicting the functional traits of individual plants [26]. Therefore, the effects of environmental factors on community leaf nutrient traits at a macro scale can reduce the impact of certain local deterministic processes (e.g., competition, mutualism).

Based on field survey data from 89 sites in China, an attempt was made to identify the key drivers influencing the leaf nutrient traits at different life forms in forest communities at the macro scale. To explain this, the following hypotheses were made: (1) there are significant differences in leaf nutrient traits at different life forms at the macro scale; and (2) climatic factors are the dominant environmental factors affecting leaf nutrient traits at different life forms of the community, while soil nutrient factors also play a coordinating role that cannot be ignored.

### **2. Results**

### *2.1. Variability in LN, LP, and N/P of Different Life Forms and Distribution Patterns in China*

Leaf nutrient traits showed significant geographical variation at different life forms, with overall higher LN and LP in northeastern China and relatively higher LP content in northwestern China (Figure 1). However, neither LN, LP, nor N/P (Figure 2) showed a significant latitudinal pattern. The LN, LP, and N/P of the herb levels showed significant differences (*p* < 0.05). The leaf nutrient traits of the arbor layer and the herb layer were significantly different (*p* < 0.001), and the herb layer had relatively higher LN and lower LP.

*Plants* **2022**, *11*, 2171 3 of 18

Leaf nutrient traits were generally significantly different between the shrub and herb levels (*p* < 0.05), but LN was not significant in the herb and shrub levels (*p* > 0.05). LP. Leaf nutrient traits were generally significantly different between the shrub and herb levels (*P* < 0.05), but LN was not significant in the herb and shrub levels (*P* > 0.05).

significant latitudinal pattern. The LN, LP, and N/P of the herb levels showed significant differences (*P* < 0.05). The leaf nutrient traits of the arbor layer and the herb layer were significantly different (*P* < 0.001), and the herb layer had relatively higher LN and lower

**Figure 1.** The distribution patterns of LN, LP, and N/P at different life forms in China with a spatial resolution of 1 × 1 km were studied by kernel density estimation. (**a**) LN of the tree levels; (**b**) LP of the tree levels; (**c**) N/P of the tree levels; (**d**) LN of the shrub levels; (**e**) LP of the shrub levels; (**f**) leaf N/P of the shrub levels; (**g**) LN of the herb levels; (**h**) herb levels LP; (**i**) herb levels N/P. **Figure 1.** The distribution patterns of LN, LP, and N/P at different life forms in China with a spatial resolution of 1 × 1 km were studied by kernel density estimation. (**a**) LN of the tree levels; (**b**) LP of the tree levels; (**c**) N/P of the tree levels; (**d**) LN of the shrub levels; (**e**) LP of the shrub levels; (**f**) leaf N/P of the shrub levels; (**g**) LN of the herb levels; (**h**) herb levels LP; (**i**) herb levels N/P.

### *2.2. Correlations between Climatic Factors and LN, LP, and N/P at Different Life Forms*

The LN and LP of the tree, shrub, or herb levels increased significantly with increasing MAT and mean annual precipitation (MAP) (Figures 3a,d and 4a,d) and decreased significantly with increasing annual sunlight duration (ASD) (Figures 3e,f and 4e). Whereas LP in both the tree and shrub levels increased significantly with increasing mean coldest monthly temperature (MCMT) and mean warmest monthly temperature (MWMT), LN in the herb level decreased significantly. The ASD and mean annual evaporation (MAE) had better predictive power for LN in the tree level compared to other climatic factors (*R <sup>2</sup>* = 0.20, *p* < 0.001; *R <sup>2</sup>* = 0.16, *p* < 0.001) and better predictive power for LP in the tree layer (*R <sup>2</sup>* = 0.46, *p* < 0.001; *R <sup>2</sup>* = 0.50, *p* < 0.001; *R <sup>2</sup>* = 0.37, *p* < 0.001) were MWMT, MAP, and ASD.

**Figure 2.** A comparison of the differences in the LN (**a**), LP (**b**), and N/P (**c**) at different life forms. (**a**) Variability of LN among the trees, shrubs, and herbs; (**b**) variability of LP among the trees, shrubs, and herbs; (**c**) variability of N/P among the trees, shrubs, and herbs. Tree represents the tree levels, Shrub represents the shrub levels, and Herb represents the herb levels. Levels are grouped where ns represents non−significant (*P* > 0.05) at the 0.05 level, \* represents *P* < 0.05, \*\* represents *P* < 0.01, \*\*\*\* represents *p* < 0.0001. **Figure 2.** A comparison of the differences in the LN (**a**), LP (**b**), and N/P (**c**) at different life forms. (**a**) Variability of LN among the trees, shrubs, and herbs; (**b**) variability of LP among the trees, shrubs, and herbs; (**c**) variability of N/P among the trees, shrubs, and herbs. Tree represents the tree levels, Shrub represents the shrub levels, and Herb represents the herb levels. Levels are grouped where ns represents non−significant (*p* > 0.05) at the 0.05 level, \* represents *p* < 0.05, \*\* represents *p* < 0.01, \*\*\*\* represents *p* < 0.0001. *Plants* **2022**, *11*, 2171 5 of 18

**Figure 3.** The general linear correlation analysis of climate factors with LN at different life forms. (**a**) General linear relationship between the MAT and LN of plants in the trees, shrubs, and herbs. (**b**) General linear relationship between the MCMT and LN of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the MWMT and LN of plants in the trees, shrubs, and herbs (**d**) General linear relationship between the MAP and LN of plants in the trees, shrubs, and herbs (**e**) General linear relationship between the ASD and LN of plants in the trees, shrubs, and herbs (**f**) General linear relationship between the MAE and LN of plants in the trees, shrubs, and herbs. MAT represents the mean annual temperature, MCMT represents the mean coldest monthly temperature, **Figure 3.** The general linear correlation analysis of climate factors with LN at different life forms. (**a**) General linear relationship between the MAT and LN of plants in the trees, shrubs, and herbs. (**b**) General linear relationship between the MCMT and LN of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the MWMT and LN of plants in the trees, shrubs, and herbs (**d**) General linear relationship between the MAP and LN of plants in the trees, shrubs, and herbs (**e**) General linear relationship between the ASD and LN of plants in the trees, shrubs, and herbs (**f**).

MWMT represents the mean warmest monthly temperature, MAP represents the mean annual precipitation, ASD represents the annual sunlight duration, and MAE represents the mean annual evaporation. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R2* represents how well the model fits the variables studied and the *P*−value rep-

resents the significance level.

General linear relationship between the MAE and LN of plants in the trees, shrubs, and herbs. MAT represents the mean annual temperature, MCMT represents the mean coldest monthly temperature, MWMT represents the mean warmest monthly temperature, MAP represents the mean annual precipitation, ASD represents the annual sunlight duration, and MAE represents the mean annual evaporation. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R 2* represents how well the model fits the variables studied and the *p*-value represents the significance level.

**Figure 4.** The general linear correlation analysis of climate factors with LP at different life forms. (**a**) General linear relationship between the MAT and LP of plants in the trees, shrubs, and herbs. (**b**) General linear relationship between the MCMT and LP of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the MWMT and LP of plants in the trees, shrubs, and herbs. (**d**) General linear relationship between the MAP and LP of plants in the trees, shrubs, and herbs. (**e**) General linear relationship between the ASD and LP of plants in the trees, shrubs, and herbs. (**f**) General linear relationship between the MAE and LP of plants in the trees, shrubs, and herbs. MAT represents the mean annual temperature, MCMT represents the mean coldest monthly temperature, MWMT represents the mean warmest monthly temperature, MAP represents the mean annual precipitation, ASD represents the annual sunlight duration, and MAE represents the mean annual evaporation. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R 2* represents how well the model fits the variables studied and the *p*-value represents the significance level.

The N/P of the herb level increased with the increase in MCMT and MWMT (Figure 5b,c), and decreased with the increase in MAT, MAP, ASD, and MAE (Figure 5a,d–f). MAT, MCMT, and MWMT (Figure 5a–c) were significantly correlated with the N/P of the shrub layer (*p* < 0.001). MAT had the best fitting effect on the N/P of the herbal layer (Figure 5b; *R <sup>2</sup>* = 0.59, *p* < 0.001).

**Figure 5.** The general linear correlation analysis of climate factors with N/P at different life forms. (**a**) General linear relationship between the MAT and N/P of plants in the trees, shrubs, and herbs. (**b**) General linear relationship between the MCMT and N/P of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the MWMT and N/P of plants in the trees, shrubs, and herbs. (**d**) General linear relationship between the MAP and N/P of plants in the trees, shrubs, and herbs. (**e**) General linear relationship between the ASD and N/P of plants in the trees, shrubs, and herbs. (**f**): General linear relationship between the MAE and N/P of plants in the trees, shrubs, and herbs. MAT represents the mean annual temperature, MCMT represents the mean coldest monthly temperature, MWMT represents the mean warmest monthly temperature, MAP represents the mean annual precipitation, ASD represents the annual sunlight duration, and MAE represents the mean annual evaporation. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R2* represents how well the model fits the variables studied and the *P*−value represents the significance level. *2.3. Effect of Soil Factors on the Relationship between LN, LP, and Leaf N/P at Different Life Forms*  **Figure 5.** The general linear correlation analysis of climate factors with N/P at different life forms. (**a**) General linear relationship between the MAT and N/P of plants in the trees, shrubs, and herbs. (**b**) General linear relationship between the MCMT and N/P of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the MWMT and N/P of plants in the trees, shrubs, and herbs. (**d**) General linear relationship between the MAP and N/P of plants in the trees, shrubs, and herbs. (**e**) General linear relationship between the ASD and N/P of plants in the trees, shrubs, and herbs. (**f**): General linear relationship between the MAE and N/P of plants in the trees, shrubs, and herbs. MAT represents the mean annual temperature, MCMT represents the mean coldest monthly temperature, MWMT represents the mean warmest monthly temperature, MAP represents the mean annual precipitation, ASD represents the annual sunlight duration, and MAE represents the mean annual evaporation. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R 2* represents how well the model fits the variables studied and the *p*-value represents the significance level.

#### The LN of the different life forms was significantly and positively correlated with the *2.3. Effect of Soil Factors on the Relationship between LN, LP, and Leaf N/P at Different Life Forms*

soil P (Figure 6b), while the trends of LP and N/P varied with the soil nutrient factors. The LN of the three life forms increased significantly with increasing soil N and P (Figure 6). Soil pH (Figure 7c) all showed a significant positive correlation (*P* < 0.001) with LP for the different life forms, whereas LN decreased significantly with increasing soil pH (Figure 6). The LP of the shrub levels increased significantly with the increasing soil N and P. The trend between the herb and shrub levels was reversed. The N/P in the herb level was significantly positively correlated with soil N (*R2* = 0.30, *P* < 0.001) (Figure 8a). Soil pH (Figure 7c) was the best predictor of LP in the herb and tree levels (*R2* = 0.25, *P* < 0.001; *R2* = 0.29, *P* < 0.001) (Figure 7c), and the best predictor of N/P in the tree level (Figure 8c) (*R2* = 0.29, *P* < 0.001) (Figure 8c). In contrast, the best prediction of N/P for the shrub level was for soil P (*R2* = 0.22, *P* < 0.001; Figure 8b). The LN of the different life forms was significantly and positively correlated with the soil P (Figure 6b), while the trends of LP and N/P varied with the soil nutrient factors. The LN of the three life forms increased significantly with increasing soil N and P (Figure 6). Soil pH (Figure 7c) all showed a significant positive correlation (*p* < 0.001) with LP for the different life forms, whereas LN decreased significantly with increasing soil pH (Figure 6). The LP of the shrub levels increased significantly with the increasing soil N and P. The trend between the herb and shrub levels was reversed. The N/P in the herb level was significantly positively correlated with soil N (*R <sup>2</sup>* = 0.30, *p* < 0.001) (Figure 8a). Soil pH (Figure 7c) was the best predictor of LP in the herb and tree levels (*R <sup>2</sup>* = 0.25, *p* < 0.001; *R <sup>2</sup>* = 0.29, *p* < 0.001) (Figure 7c), and the best predictor of N/P in the tree level (Figure 8c) (*R <sup>2</sup>* = 0.29, *p* < 0.001) (Figure 8c). In contrast, the best prediction of N/P for the shrub level was for soil P (*R <sup>2</sup>* = 0.22, *p* < 0.001; Figure 8b).

*Plants* **2022**, *11*, 2171 8 of 18

**Figure 6.** The general linear analysis of soil factors with different life forms of LN. (**a**) General linear relationship between the soil N and LN of plants in the trees, shrubs, and herbs. (**b**) General linear relationship between the soil P and LN of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the soil pH and LN of plants in the trees, shrubs, and herbs. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R2* represents how well the model fits the variables studied and *P*−value represents the level of significance. **Figure 6.** The general linear analysis of soil factors with different life forms of LN. (**a**) General linear relationship between the soil N and LN of plants in the trees, shrubs, and herbs. (**b**) General linear relationship between the soil P and LN of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the soil pH and LN of plants in the trees, shrubs, and herbs. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R 2* represents how well the model fits the variables studied and *p*-value represents the level of significance. **Figure 6.** The general linear analysis of soil factors with different life forms of LN. (**a**) General linear relationship between the soil N and LN of plants in the trees, shrubs, and herbs. (**b**) General linear relationship between the soil P and LN of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the soil pH and LN of plants in the trees, shrubs, and herbs. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R2* represents how well the model fits the variables studied and *P*−value represents the level of significance.

relationship between the soil P and LP of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the soil pH and LP of plants in the trees, shrubs, and herbs. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R2* represents how well the model fits the variables studied and the *P*−value represents the level of significance. **Figure 7.** The general linear analysis of soil factors with different life forms of LP. (**a**) General linear relationship between the soil N and LP of plants in the trees, shrubs, and herbs (**b**) General linear relationship between the soil P and LP of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the soil pH and LP of plants in the trees, shrubs, and herbs. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R2* represents how well the model fits the variables studied and the *P*−value represents the level of significance. **Figure 7.** The general linear analysis of soil factors with different life forms of LP. (**a**) General linear relationship between the soil N and LP of plants in the trees, shrubs, and herbs (**b**) General linear relationship between the soil P and LP of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the soil pH and LP of plants in the trees, shrubs, and herbs. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R 2* represents how well the model fits the variables studied and the *p*-value represents the level of significance. *Plants* **2022**, *11*, 2171 9 of 18

**Figure 8.** The general linear analysis of soil factors with different life forms of N/P. (**a**) General linear relationship between the soil N and N/P of plants in the trees, shrubs, and herbs (**b**) General linear **Figure 8.** The general linear analysis of soil factors with different life forms of N/P. (**a**) General linear

relationship between the soil P and N/P of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the soil pH and N/P of plants in the trees, shrubs, and herbs. The red line

*2.4. Climatic and Soil Factors Dominate Changes in the Functional Traits of Different Commu-*

de = 42.2%, 41%, 22.8%) and LP (Figure 10d−f; de = 17.2%, 17.1%, 50.6%).

We analyzed the effects of environmental factors on the leaf LN, LP, and N/P at different life forms based on a generalized additive model with non−metric multidimensional scaling (NMDS) ranking. Overall, the leaf N/P showed strong environmental plasticity, and environmental factors generally explained the spatial variation in leaf N/P more than LN and LP at different life forms (Figures 9–11). Climate factors played a greater role than the soil nutrient factors in shaping the spatial variation in the leaf LN (Figure 9d−f;

significance.

*nities* 

relationship between the soil N and N/P of plants in the trees, shrubs, and herbs (**b**) General linear relationship between the soil P and N/P of plants in the trees, shrubs, and herbs. (**c**) General linear relationship between the soil pH and N/P of plants in the trees, shrubs, and herbs. The red line represents the tree level, the green line is the shrub level, and the blue line is the herb level. *R 2* represents how well the model fits the variables studied and the *p*-value represents the level of significance.

### *2.4. Climatic and Soil Factors Dominate Changes in the Functional Traits of Different Communities*

We analyzed the effects of environmental factors on the leaf LN, LP, and N/P at different life forms based on a generalized additive model with non−metric multidimensional scaling (NMDS) ranking. Overall, the leaf N/P showed strong environmental plasticity, and environmental factors generally explained the spatial variation in leaf N/P more than LN and LP at different life forms (Figures 9–11). Climate factors played a greater role than the soil nutrient factors in shaping the spatial variation in the leaf LN (Figure 9d−f; de = 42.2%, 41%, 22.8%) and LP (Figure 10d−f; de = 17.2%, 17.1%, 50.6%). *Plants* **2022**, *11*, 2171 10 of 18

**Figure 9.** The NMDS ranking of climatic and soil factors with different life forms of LN. Value of de represents the deviation explained by the corresponding model. (**a**) NMDS ranking of soil factors **Figure 9.** The NMDS ranking of climatic and soil factors with different life forms of LN. Value of de

with tree levels LN; (**b**) NMDS ranking of soil factors with shrub levels LN; (**c**) NMDS ranking of soil factors with herb levels LN; (**d**) NMDS ranking of climatic factors with tree levels LN; (**e**) NMDS

ranking of the sum of soil factors and climate factors and shrub levels LN; (**i**) NMDS ranking of the sum of soil factors and climate factors and herb levels LN. Trait stacking indicates that abiotic factors, indicated by points on the NMD, are associated with higher or lower trait values, consistent with a colored trait gradient. Note that if the relationship between the LN and abiotic factors is linear, the gradient splines will be parallel. Nonlinear relationships between LN and abiotic factors are repre-

sented by curve splines.

represents the deviation explained by the corresponding model. (**a**) NMDS ranking of soil factors with tree levels LN; (**b**) NMDS ranking of soil factors with shrub levels LN; (**c**) NMDS ranking of soil factors with herb levels LN; (**d**) NMDS ranking of climatic factors with tree levels LN; (**e**) NMDS ranking of climatic factors with shrub levels LN; (**f**) NMDS ranking of climate factors and herb levels LN; (**g**) NMDS ranking of the sum of soil factors and climate factors and tree levels LN; (**h**) NMDS ranking of the sum of soil factors and climate factors and shrub levels LN; (**i**) NMDS ranking of the sum of soil factors and climate factors and herb levels LN. Trait stacking indicates that abiotic factors, indicated by points on the NMD, are associated with higher or lower trait values, consistent with a colored trait gradient. Note that if the relationship between the LN and abiotic factors is linear, the gradient splines will be parallel. Nonlinear relationships between LN and abiotic factors are represented by curve splines. *Plants* **2022**, *11*, 2171 11 of 18

**Figure 10.** The NMDS ranking of climatic and soil factors with different life forms of LP. Value of de represents the deviation explained by the corresponding model. (**a**) NMDS ranking of soil factors with tree levels LP; (**b**) NMDS ranking of soil factors with shrub levels LP; (**c**) NMDS ranking of soil factors with herb levels LP; (**d**) NMDS ranking of climatic factors with tree levels LP; (**e**) NMDS ranking of climatic factors with shrub levels LP; (**f**) NMDS ranking of climate factors and herb levels LP; (**g**) NMDS ranking of the sum of soil factors and climate factors and tree levels LP; (**h**): NMDS ranking of the sum of soil factors and climate factors and shrub levels LP; (**i**): NMDS ranking of the **Figure 10.** The NMDS ranking of climatic and soil factors with different life forms of LP. Value of de represents the deviation explained by the corresponding model. (**a**) NMDS ranking of soil factors with tree levels LP; (**b**) NMDS ranking of soil factors with shrub levels LP; (**c**) NMDS ranking of soil factors with herb levels LP; (**d**) NMDS ranking of climatic factors with tree levels LP; (**e**) NMDS ranking of climatic factors with shrub levels LP; (**f**) NMDS ranking of climate factors and herb levels

sum of soil factors and climate factors and herb levels LP. Trait stacking indicates that abiotic factors, indicated by points on the NMD, are associated with higher or lower trait values, consistent with a

by curve splines.

LP; (**g**) NMDS ranking of the sum of soil factors and climate factors and tree levels LP; (**h**): NMDS ranking of the sum of soil factors and climate factors and shrub levels LP; (**i**): NMDS ranking of the sum of soil factors and climate factors and herb levels LP. Trait stacking indicates that abiotic factors, indicated by points on the NMD, are associated with higher or lower trait values, consistent with a colored trait gradient. Note that if the relationship between LP and abiotic factors is linear, the gradient splines will be parallel. Nonlinear relationships between LP and abiotic factors are represented by curve splines. *Plants* **2022**, *11*, 2171 12 of 18

**Figure 11.** The NMDS ranking of climatic and soil factors with different life forms of N/P. Value of de represents the deviation explained by the corresponding model. (**a**) NMDS ranking of soil factors with tree levels N/P; (**b**) NMDS ranking of soil factors with shrub levels N/P; (**c**) NMDS ranking of soil factors with herb levels N/P; (**d**) NMDS ranking of climatic factors with tree levels N/P; (**e**) NMDS ranking of climatic factors with shrub levels N/P; (**f**) NMDS ranking of climate factors and herb levels N/P; (**g**) NMDS ranking of the sum of soil factors and climate factors and tree levels N/P; (**h**) NMDS ranking of the sum of soil factors and climate factors and shrub levels N/P; (**i**) NMDS ranking of the sum of soil factors and climate factors and herb levels N/P. Trait stacking indicates that abiotic factors, indicated by points on the NMD, are associated with higher or lower trait values, consistent with a colored trait gradient. Note that if the relationship between N/P and abiotic factors is linear, the gradient splines will be parallel. Nonlinear relationships between N/P and abiotic factors are represented by curve splines. **Figure 11.** The NMDS ranking of climatic and soil factors with different life forms of N/P. Value of de represents the deviation explained by the corresponding model. (**a**) NMDS ranking of soil factors with tree levels N/P; (**b**) NMDS ranking of soil factors with shrub levels N/P; (**c**) NMDS ranking of soil factors with herb levels N/P; (**d**) NMDS ranking of climatic factors with tree levels N/P; (**e**) NMDS ranking of climatic factors with shrub levels N/P; (**f**) NMDS ranking of climate factors and herb levels N/P; (**g**) NMDS ranking of the sum of soil factors and climate factors and tree levels N/P; (**h**) NMDS ranking of the sum of soil factors and climate factors and shrub levels N/P; (**i**) NMDS ranking of the sum of soil factors and climate factors and herb levels N/P. Trait stacking indicates that abiotic factors, indicated by points on the NMD, are associated with higher or lower trait values,

**3. Discussion** 

consistent with a colored trait gradient. Note that if the relationship between N/P and abiotic factors is linear, the gradient splines will be parallel. Nonlinear relationships between N/P and abiotic factors are represented by curve splines.

### **3. Discussion**

The leaf nutrient traits (LN, LP, and N/P) of different plant forms showed significant geographical differences [27], and some studies have found that the LN of shrub levels is significantly higher than that of the tree and herb levels in northwestern China, while the LP of shrub levels is significantly lower than that of herb levels [28], and that the LP and N/P of shrubs are higher than those of herbs in the desertification zone of southwest Hunan [29]. Within a given community, differences in the ecological niches of species and differences in resource use patterns can lead to differences in the leaf nutrient traits between species [17]. Leaf nutrient traits are widely used to quantify plant adaptations to the environment, and even to predict ecosystem function [6,30]. Thus, the differences in LN and LP of different life forms of plants result from differences in plant ecological niches. Compared to herb plants, woody plants have lower LP, which is aligned with their low growth rate strategy [31]. Quantifying the geographical distribution patterns of leaf nutrient traits is important for our scientific assessment of forest development dynamics [24].

Temperature can have a significant effect on the LN, LP, and N/P [32], and the temperature−biogeochemical hypothesis also suggests that LN and LP increase monotonically with temperature on a global scale [33]. It has been found that LN and LP are strongly correlated with effective soil N and P content [32], and that climatic factors can act directly on soil nutrients [25]. Soil enzymes are beneficial to maintain soil fertility, and increasing the soil temperature and moisture can promote the activity of soil enzymes and significantly improve the decomposition efficiency of soil organic matter [34]. Low temperatures not only have an inhibitory effect on organic matter decomposition and mineralization, but also limit microbial activity and affect the decomposition of apoplastic matter, thus reducing the availability of soil N and P [33]. Conversely, as the temperature increases, soil microbial respiration is enhanced, the efficiency of organic matter mineralization and decomposition increases, and the effective soil N, and P content rises [35]. Therefore, MAT is a key limiting factor for changes in leaf nutrient traits [12,13]. It has also been found that ASD significantly affects LN and LP [36]. Plant N and P elements are involved in the light reaction process and have an impact on photosynthesis, which is enhanced with increasing light duration and LN and LP are heavily utilized [37]. The nutrient content of plant leaves gradually decreased with increasing light time. Our results demonstrate the important influence of temperature on the nutrients of different life forms of leaves.

Precipitation is also one of the main factors of nutrient traits in plant leaves [14]. Previous studies have found that in water−limited ecosystems, plant N and P uptake of elements is mainly limited by water [38]. On one hand, most soil water comes from rainfall, and soil water content can directly affect the activity of soil enzymes, thereby affecting the release of soil nutrients [34]. On the other hand, rainfall promotes plant N uptake by promoting soil element availability, and plant N uptake is positively correlated with soil N availability [39]. In addition, rainfall had a significant effect on the soil microbial respiration and accelerated the leaching and transformation of P, which contributed to the transformation of soil P and thus increased the effective soil P content. With increased precipitation, the effectiveness of the soil nutrient factors increased, and plants were able to take up more N and P, significantly increasing the LN and LP [40]. Thus, precipitation is an important effect on the spatial and temporal distribution patterns of plant nutrient traits.

Soil N and P are basic nutrients for plant growth and are closely related to plant leaf nutrient, and most of the LN and LP of plants originate from soil [9]. The availability of soil N and P elements directly determines the growth and development process of plants, and available soil nutrients are positively correlated with leaf nutrients [18]. It was found that fertilizing the soil with N and P was beneficial to maintaining and improving the activity of soil enzymes and releasing more micronutrients and nutrients [41]. The higher the effective soil N and P content, the higher the plant root N and P uptake efficiency, and the higher the LN and N/P [22]. Soil microbial activity is positively correlated with soil pH, especially the activity of soil microbial enzymes related to the breakdown of soil nutrients. As the pH increases, soil microbial metabolic activity decreases, inhibiting the soil N mineralization and plant uptake of soil N, and LN is subsequently reduced [42]. However, it was also found that soil microbial biomass increased with the increasing soil pH, effectively driving organic P mineralization and the dissolution of fixed P [43], and LP in the arboreal and herbaceous layers increased. Thus, soil nutrient factors significantly influenced the changes in the leaf nutrient traits.

At the macro scale, climate factors have a stronger ability to shape leaf traits than soil factors [44]. Studies have shown that plant leaves are more sensitive to changes in climate and that plants can respond to changing climate by changing their traits [1]. Soil nutrient elements are closely related to leaf nutrient elements and provide many nutrients for plant growth and development, but the availability of soil nutrients is limited by temperature and precipitation [33]. Thus, climatic factors play a stronger role in shaping leaf nutrient distribution patterns than soil factors [45], but the influence of soil nutrient factors on leaf traits is not negligible.

The research area of this thesis included most of areas, which together with the relatively large amount of data ensures the plausibility of our experimental results. In addition, the results of this study are presented using a variety of data analysis methods. Overall, our results quantify the relative contribution of different environmental factors in shaping the functional traits of leaf nutrients and are important for a better understanding of the impact of global climate change on plant physiology [10]. In future studies, more biological factors such as community biodiversity and stand density need to be further explored for their effects on the distribution patterns of community leaf nutrients.

### **4. Materials and Methods**

### *4.1. Study Area and Sample Data*

The relatively wide distribution of forests in China facilitates studies on large scales. Data from 818 forest at 89 sites plots surveyed between 2005 and 2020 were used to investigate the spatial distribution and driving factors of LN, LP, and N/P at different life forms (tree levers, shrub levels, herb levels) of forest communities (Figure S1). The study sites ranged from 19.1◦ to 53.5◦ N and 79.7◦ to 129.3◦ E.

### *4.2. Functional Data*

Community nutrient profiles can better reflect the adaptability of local vegetation to different environments than individual nutrient profiles [46]. Data for the selected nutrient profiles in this study include LN, LP, and N/P.

In this experiment, LN was determined using the national standard method, the Kjeldahl method [47], LP was determined using the vanadium−molybdenum yellow absorbance method [48], and N/P was equal to LN/LP.

Theoretical risk in predicting functional traits in individual plants is due to intraspecific and interspecific struggles [49]. Community−weighted mean traits (*CWM*<sup>i</sup> ) represents the forest mean trait values.

$$\text{CWM}\_{\text{i}} = \frac{\sum\_{\text{i}}^{\text{n}} D\_{\text{i}} \times \text{Trait}}{\sum\_{\text{i}}^{\text{n}} D\_{\text{i}}}$$

where *CWM*<sup>i</sup> represents the community−weighted functional trait identity value and *D*<sup>i</sup> represents the abundance of dominant tree species. *Trait*<sup>i</sup> represents the selected functional trait [50].

### *4.3. Environmental Data*

A study found that changes in the functional traits of building blocks were related to climate change [51]. The mean annual temperature (MAT), mean coldest monthly temperature (MCMT), mean warmest monthly temperature (MWMT), and mean annual precipitation (MAP) were extracted from the WorldClim global climate layer at a spatial resolution of 1 km. The mean annual evaporation (MAE) was taken from the Climate Data Center of the China Meteorological Administration (http://data.cma.cn/site/index. html) (accessed on 1 April 2021), sunlight. The annual sunlight duration (ASD), a key factor in photosynthesis [52], is from the Climate Data Center of the China Meteorological Administration (http://data.cma.cn/site/index.html) (accessed on 1 April 2021). Soil pH, soil N (http://www.csdn.store) (accessed on 1 April 2021), and soil P (https://www.osgeo. cn/data/wc137) (accessed on 1 April 2021) in the top 30 cm of soil were extracted from a 250 m resolution grid.

### *4.4. Data Analysis*

Data analysis of the LN, LP, and N/P data from our study conformed to a normal distribution (Figure S2). A significant difference test at the 0.05 significance level was used to test for significant differences in LN, LP, and N/P between the tree, shrub, and herb levels (Figure 2). Significant differences were analyzed using the R package agricolae (version 4.1.0, R Core Team, 2020).

The extent to which environmental factors explain LN, LP, and N/P was investigated using a linear regression model in the R package agricolae (version 4.1.0, R Core Team, 2020), where *R 2* represents how well the model fits the variables studied.

The generalized additive model (GAM) is used to test the effects of various environmental factors on the functional traits of leaves, with data for the LN, LP, and N/P consisting of parametric and non−parametric components to reduce the model risk associated with linear models [49]. First, the GAM method was used to select the key influencing factors based on a significant difference test at the 0.05 level of significance. Then, a GAM model was constructed to measure the relationship between the environmental factors and LN, LP, and N/P and used non−metric multidimensional scaling analysis (NMDS) to reflect the results for GAM [53].

$$\lg[\mathcal{E}(\mathbf{Y}|\mathbf{X})] = \sum\_{\mathbf{i}} \mathfrak{H}\mathbf{i}\mathbf{X}\mathbf{i} + \sum\_{\mathbf{j}} \mathrm{fi}(\mathbf{X}\mathbf{i}) + \varepsilon\_{\mathbf{j}}$$

where g(•) represents the connection function, the form depends on the specific form, which can be interpreted as the Y−variable distribution. Є is the random error term, which can be interpreted as the variable connection with the normal distribution function name identity, and the connection function has the form g(u) = u, u = E(Y|X), E(є|X) = 0. X<sup>i</sup> is the explanatory variable that strictly follows the parametric form, β<sup>i</sup> is the corresponding parameter, and f<sup>j</sup> (X<sup>j</sup> ) is the corresponding explanatory variable that follows the nonparametric form of the smoothing function. In our study, the spline smoothing function S(•) was selected for fitting, thin−plate spline smoothing was selected for function fitting between different nodes, and each smoothing function S(•) was estimated using penalized least squares [50].

### **5. Conclusions**

This study used data at 89 sites from 818 forest plots across China over 15 years from 2005 to 2020 to verify the relative roles of climatic and soil factors in shaping LN, LP, and N/P at different life forms. Our results confirm that LN, LP, and N/P differ significantly between life forms and that climatic and soil factors in the community habitat jointly influence the distribution patterns of the N/P nutrient profiles, with MAT being the dominant factor influencing LN, LP, and N/P at different life forms. Climatic factors also indirectly influence the leaf nutrient traits by affecting soil nutrients, but climatic factors are more influential than soil factors in shaping the distribution patterns of LN, LP, and N/P. This study has important implications for our understanding of the distribution patterns of the plant N and P nutrient profiles at different life forms in the context of climate change and ecosystem modeling.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11162171/s1, Figure S1: Sample site map of the distribution of the sample sites used for this experimental study in China. The sample sites ranged from 19.1◦ N to 53.5◦ N and 79.7◦ E to 129.3◦ E.; Figure S2: The normal distribution charts for the LN, LP, and N/P data.

**Author Contributions:** X.W. and J.W. contributed equally to this work. Conceptualization, J.G.; Experiment implementation, J.G.; Validation, J.G.; Formal analysis, X.W. and J.W.; Writing—original draft preparation, X.W.; Writing—review and editing, X.W., J.W., C.L., L.L. and L.Z.; Visualization, J.G.; Project administration, J.G. and H.Z.; Funding acquisition, J.G. and H.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** The study was supported by the Xinjiang Normal University Landmark Achievements Cultivation Project, China (grant number: no number), the Scientific Research Program of Colleges and Universities in Xinjiang (No. XJEDU2021I023), the Natural Science Foundation of China (No. 32160074), and the Open Project of Key Laboratory in Xinjiang (No. 2020D4010).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

### **References**


**Lin Wang <sup>1</sup> , Qinyue Miao <sup>1</sup> , Ülo Niinemets 2,3,\*, Johan Gielis <sup>4</sup> and Peijian Shi 1,\***


**Abstract:** Many geometries of plant organs can be described by the Gielis equation, a polar coordinate

equation extended from the superellipse equation, *r* = *a* h cos *<sup>m</sup>* <sup>4</sup> ϕ *<sup>n</sup>*<sup>2</sup> + 1 *k* sin *<sup>m</sup>* <sup>4</sup> ϕ *n*3 i−1/*n*<sup>1</sup> . Here, *r* is the polar radius corresponding to the polar angle ϕ; *m* is a positive integer that determines the number of angles of the Gielis curve when ϕ ∈ [0 to 2π); and the rest of the symbols are parameters to be estimated. The pentagonal radial symmetry of calyxes and corolla tubes in top view is a common feature in the flowers of many eudicots. However, prior studies have not tested whether the Gielis equation can depict the shapes of corolla tubes. We sampled randomly 366 flowers of *Vinca major* L., among which 360 had five petals and pentagonal corolla tubes, and six had four petals and quadrangular corolla tubes. We extracted the planar coordinates of the outer rims of corolla tubes (in top view) (ORCTs), and then fitted the data with two simplified versions of the Gielis equation with *k* = 1 and *m* = 5: *r* = *a* h cos 5 <sup>4</sup>ϕ *n*2 + sin 5 <sup>4</sup>ϕ *n*3 i−1/*n*<sup>1</sup> (Model 1), and *r* = *a* h cos 5 <sup>4</sup>ϕ *n*2 + sin 5 <sup>4</sup>ϕ *n*2 i−1/*n*<sup>1</sup> (Model 2). The adjusted root mean square error (RMSEadj) was used to evaluate the goodness of fit of each model. In addition, to test whether ORCTs are radially symmetrical, we correlated the estimates of *n*<sup>2</sup> and *n*<sup>3</sup> in Model 1 on a log-log scale. The results validated the two simplified Gielis equations. The RMSEadj values for all corolla tubes were smaller than 0.05 for both models. The numerical values of *n*<sup>2</sup> and *n*<sup>3</sup> were demonstrated to be statistically equal based on the regression analysis, which suggested that the ORCTs of *V. major* are radially symmetrical. It suggests that Model 1 can be replaced by the simpler Model 2 for fitting the ORCT in this species. This work indicates that the pentagonal or quadrangular corolla tubes (in top view) can both be modeled by the Gielis equation and demonstrates that the pentagonal or quadrangular corolla tubes of plants tend to form radial symmetrical geometries during their development and growth.

**Keywords:** flower geometry; Gielis equation; model complexity; natural geometries; planar coordinates; polygonal structure; radial symmetry

### **1. Introduction**

In geometry, the circle and ellipse have been demonstrated to be expressed as two special cases of the superellipse [1]. As a direct extension of the superellipse equation in the polar coordinate system, Gielis [2,3] created a highly versatile equation to reflect natural geometries, especially symmetrical geometries. We refer to it as the Gielis equation hereinafter. In the past decade, many studies have been carried out to examine the validity of the Gielis equation in fitting actual biological geometries. A simplified version of this

**Citation:** Wang, L.; Miao, Q.; Niinemets, Ü.; Gielis, J.; Shi, P. Quantifying the Variation in the Geometries of the Outer Rims of Corolla Tubes of *Vinca major* L. *Plants* **2022**, *11*, 1987. https://doi.org/ 10.3390/plants11151987

Academic Editor: Yasutomo Hoshika

Received: 4 July 2022 Accepted: 26 July 2022 Published: 30 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

equation was used by Shi et al. [4,5] and Lin et al. [6] to fit the boundary data of crosssections of tree rings for five species of conifers and bamboo leaves for 46 bamboo species. These studies verified the potential of the Gielis equation to describe the shapes of tree rings and bamboo leaves and demonstrated that the Gielis equation is a valid scientific method for quantitative characterization of the size and shape of widely differing planar biological objects.

Shi et al. [7] proposed a twin version of the Gielis equation by introducing a link function and found that the twin Gielis equation was superior in depicting the shapes of some sea stars. Tian et al. [8] used the Gielis equation to fit the seed projections (in side view) of two *Gingko biloba* cultivars, and used it to quantify the morphological differences between the two cultivars. Li et al. [9] compared the original Gielis equation with its twin version in describing the planar projections of *Koelreuteria paniculata* fruits (in top view) and demonstrated that the two versions of the Gielis equation both can model the shapes of the vertical fruit projections well. Nevertheless, the twin Gielis equation predicted an axial symmetry of *K. paniculata* fruits, while the original Gielis equation tends to overfit the data. In addition, a recent study shows that a simplified Gielis equation can describe all existing egg shapes of birds and has a better goodness of fit than other egg shape models [10]. Flowers can also be modelled in the same way [2,3,11–13], but the capacity of the Gielis equation to simulate flower shape has not been quantitatively studied.

Here, we focused on the geometries of the outer rims of corolla tubes of *Vinca* flowers that are representative to the five-petal flowers with a fused pentagonal base (corolla top). The genus *Vinca* is native to western Mediterranean Europe, Asia Minor and Northern Africa, but it has been introduced as ornamental to all continents and has naturalized in many sites. It is a small evergreen ground cover plant and in addition to sexual reproduction, spreads via stolons. Taxonomically, *Vinca* belongs to the Apocynaceae family, one of the five families within the order Gentianales; together with the orders Solanales and Lamiales they form the Lamiids [14]. Lamiids are characterized by late sympetaly, the fusion of stamen filaments with the corolla tube and opposite leaves [15]. The Apocynaceae have a conserved architecture of highly synorganized flowers, and within this family *Vinca* L. is the type genus of the tribe Vinceae, in particular, of the subtribe Vincinae. The corolla is infundibuliform, and the lobe aestivation is sinistrorse [15].

The petals of *Vinca* are fused at their bottom, forming a corolla tube. When the flower opens, the distal, unfused parts of the petals fold back in a plane, whereas the upper part of the corolla tube formed by the fused parts of the petals becomes clearly delineated. The upper ridge or rim of the corolla tube has a clear pentagonal symmetry (although quadrangular or hexagonal symmetry may occur). In contrast to the purple color of the free petals and the base of the tube, the upper rim of the floral tube is white. The purple petals exhibit a high ultraviolet reflectance, whereas the corolla tube, in particular, the upper rim of the flower, strongly absorbs ultraviolet (Ultraviolet Flowers: Vinca minor. Available online: http://www.naturfotograf.com/UV\_VINC\_MIN.htm (accessed on 1 July 2022)).

Most flowers in this species have five petals, but flowers with four or six petals occur rarely (Figure 1). Five petals correspond to a corolla tube of pentagonal symmetry and a pentagonal rim, and four petals correspond to a quadrangular rim (Figure 1). In the remainder of the paper, we use the term the outer rim of the corolla tube (in top view) (ORCT) to denote this polygonal structure and characterize the geometry of the outer rim of the corolla tube in top view (i.e., represented by a vertical projection). The two types of ORCTs (i.e., pentagonal and quadrangular rims) both seem to exhibit a radial symmetry, whereas the distal ends of the petals are rotated counterclockwise, relative to the corolla tube. The main reason is the sinistrorse aestivation of the flowers. The *Vinca* flowers display contorted aestivation of the corolla, and each petal is asymmetric. However, the entire corolla exhibits a rotational symmetry (Figure 1).

**Figure 1.** Representative flowers of *Vinca major* L. with five (**a**) and with four petals (**b**).

Pentagonal symmetry is a general condition in Eudicots [16], but a clear pentagonal shape, as in the corolla tube of *Vinca*, is uncommon. In trumpet, campanulate or salverform corolla tubes, the transition from fused proximal parts of the petals to the free parts is gradually curved. In *Vinca*, on the other hand, the plane formed by the free ends of the petals is almost perpendicular to the corolla tube. Although the ORCT is never completely flat, a projection using a photography or an image scanner can be obtained, making a 2D quantitative study possible. In this study, 360 flowers with five petals and six flowers with four petals from *V. major* were sampled to examine whether the ORCTs follow the Gielis equation, and to evaluate whether a deviation from pure rotational symmetry can be found.

### **2. Materials and Methods**

### *2.1. Flower Sampling and Image Processing*

*Vinca major* flowers were randomly sampled at the Nanjing Forestry University campus (118◦4803500 E, 32◦0406700 N), Nanjing, China from 7–23 April 2022 when the peak blooming occurred. To keep the flowers fresh, each sampled flower was placed in a 10 mL beaker with 1–2 mL water until its image was scanned. The flowers were scanned by an Epson photo scanner (V550, Epson, Batam, Indonesia) at a resolution of 2400 dpi, then their images were converted into black-white images after being cropped and saved as bmp (Bitmap) format using Adobe Photoshop CS2 (version 9.0; Adobe, San Jose, CA, USA; http://www.adobe.com/products/photoshop.html, accessed on 1 July 2022).

### *2.2. Data Acquisition*

To extract the planar coordinates of the outer rim of the corolla tube (in top view) (ORCT), we used MATLAB (version ≥ 2009a; MathWorks, Natick, MA, USA) with a program developed by refs. [4,17,18]. We fitted the coordinate data of each ORCT using the 'biogeom' package (version 1.0.5) [19] based on R (version 4.2.0) [20].

### *2.3. Models*

Gielis [2] proposed a polar coordinate equation to describe natural shapes:

$$r(\varphi) = \left( \left| \frac{1}{A} \cos \left( \frac{m}{4} \varphi \right) \right|^{n\_2} + \left| \frac{1}{B} \sin \left( \frac{m}{4} \varphi \right) \right|^{n\_3} \right)^{-\frac{1}{n\_1}} \tag{1}$$

where *r* and ϕ are the polar radius and polar angle, respectively; *A*, *B*, *n*1, *n*<sup>2</sup> and *n*<sup>3</sup> are parameters to be fitted; and *m* is a positive integer that determines the number of angles of the Gielis curve within [0, 2π). Given that the ORCT can be pentagonal or quadrangular, *m* was either five or four in this study.

Equation (1) can be rewritten as [7,8]:

$$r(\varphi) = a \left( \left| \cos \left( \frac{m}{4} \varphi \right) \right|^{n\_2} + \left| \frac{1}{k} \sin \left( \frac{m}{4} \varphi \right) \right|^{n\_3} \right)^{-\frac{1}{n\_1}} \tag{2}$$

where *a* = *A <sup>n</sup>*2/*n*<sup>1</sup> and *k* = *B*/*A <sup>n</sup>*2/*n*<sup>3</sup> . In this work, *k* was set to 1 because the ORCT of *V. major* exhibits radial symmetry. Most flowers have five petals, so *m* was set to 5, resulting in:

$$r(\varphi) = a \left( \left| \cos \frac{5}{4} \varphi \right|^{n\_2} + \left| \sin \frac{5}{4} \varphi \right|^{n\_3} \right)^{-\frac{1}{n\_1}} \tag{3}$$

To test whether a radial symmetrical version of the Gielis equation is applicable to the ORCT, we set *n*<sup>2</sup> = *n*<sup>3</sup> in Equation (3) so that we have:

$$r(\varphi) = a \left( \left| \cos \frac{5}{4} \varphi \right|^{n\_2} + \left| \sin \frac{5}{4} \varphi \right|^{n\_2} \right)^{-\frac{1}{n\_1}} \tag{4}$$

We refer to Equations (3) and (4) as Models 1 and 2 for convenience hereinafter.

### *2.4. Model Fitting and Data Analysis*

We used the Nelder–Mead optimization method [21] to minimize the residual sum of square (RSS) between the observed and predicted polar radii, and obtained the estimated values of the parameters in Equations (3) and (4):

$$\text{RSS} = \sum\_{i=1}^{N} \left( r\_i - \mathfrak{r}\_i \right)^2 \tag{5}$$

where *r<sup>i</sup>* represents the observed distance from the polar point to the *i*-th point on a scanned ORCT; *r*ˆ*<sup>i</sup>* represents the predicted distance from the polar point to the *i*-th point on the predicted ORCT based on Equation (3) or Equation (4); and *N* represents the number of data points on a scanned ORCT.

Additionally, the root-mean-square error (RMSE) was calculated to reflect the goodness of fit:

$$\text{RMSE} = \sqrt{\text{RSS}/N} \tag{6}$$

However, given the influence of the ORCT size (area) on absolute values of RMSE, we used the adjusted RMSE (RMSEadj) [7,9,22]. In that case, we can directly compare the differences in the model goodness of fits among different ORCT sizes.

$$\text{RMSE}\_{\text{adj}} = \sqrt{\frac{\text{RSS}/N}{\text{S}/\pi}}\tag{7}$$

where *S* represents the area of an ORCT. The RMSEadj represents the ratio of the mean absolute deviation (between the observed and predicted radii from the polar point to the ORCT) to the radius of a hypothetical circle whose area equals to that of the ORCT projection, which can standardize the prediction error regardless of the ORCT size.

There are still four parameters in Model 1, and the question is whether a less complex model can produce an approximate goodness of fit with fewer parameters. On the one hand, the difference between *n*<sup>2</sup> and *n*<sup>3</sup> in Model 1 determines the extent of symmetry. If the difference is equal to zero (i.e., Model 2), it can produce a perfectly axial symmetrical curve for the ORCT; the larger the difference, the worse the extent of the rim's symmetry. On the other hand, the ORCTs of interest appear visually to be axial symmetrical pentagons, so it is necessary to test whether Model 1 (with four model parameters) can be simplified to Model 2 (with three model parameters). Therefore, we performed a linear regression on the estimated values of *n*<sup>3</sup> and *n*<sup>2</sup> for all of the samples. To stabilize the variance of the estimated values of the two parameters and to normalize the data, the log-transformation was used [23,24]:

$$y = \alpha + \beta \ge \tag{8}$$

where *y* = ln *n*ˆ <sup>3</sup> and *x* = ln *n*ˆ <sup>2</sup>, where the circumflex represents the estimated value. The parameters α and β were estimated using reduced major axis regression protocols [25,26]. The bootstrap percentile method [27,28] was used to calculate 95% confidence intervals (CIs) of the intercept and the slope of the regression line. If the CI of the intercept includes zero and the CI of the slope includes one, it can indicate that *n*<sup>3</sup> is not significantly different from *n*2, which means that Model 2 is superior to Model 1 due to reduced complexity. It also suggests that ORCTs tend to be perfectly axially symmetrical. If *n*<sup>3</sup> is not statistically significant from *n*2, it is unnecessary to compare the RMSEadj values of the two models. It is apparent that an additional parameter can increase the goodness of fit, but it might result from the overfitting to the data with a certain measurement error. In other words, the flexibility in curve fitting of Model 1 may cause an incorrect parameter estimation due to the measurement errors. As a rule of thumb, a ≤ 0.05 RMSEadj can reflect the validity of a model in curve fitting.

To compare the goodness of fit between Models 1 and 2, a paired sample *t*-test was used to compare the average values of the RMSEadj values of Models 1 and 2 at 0.05 significance level. The log-transformed values of adjusted RMSEs were used to normalize the data of adjusted RMSEs.

In addition, we used nonlinear least-squares to fit the ORCTs by minimizing the sum of squares (RSS) between the observed and predicted polar radii. The distribution of the residuals in nonlinear least-squares is usually hypothesized to be normal, but the hypothesis in nonlinear least-squares seldom holds true. In that case, we calculated the mean and corresponding 95% CI of residuals to examine whether the mean equals 0 and the corresponding 95% CI includes 0. We also calculated the skewness (*S<sup>k</sup>* ) of residuals to see whether the skewness seriously deviates from zero. For a normal distribution, the skewness equals zero. A positive skewness represents a right skewed distribution curve, and a negative skewness represents a right-skewed distribution curve [17]:

$$S\_k = E\left[\left(\frac{z-\mu}{\sigma}\right)^3\right] \tag{9}$$

where *z* represents the residual between the *i*-th observed and predicted radii; µ represents the mean of residuals; and σ represents the standard error of residuals. We used the Shapiro–Wilk test [29] to test the normality of residuals only as a reference.

### **3. Results**

Both Models 1 and 2 provided generally a good representation of the ORCT (Figure 2 for sample fits). The RMSEadj values of pentagonal ORCTs predicted by Models 1 and 2 ranged from 0.0116 to 0.0392 and from 0.0118 to 0.0481, respectively (Figure 3a). This verified the validities of the two models in describing the shapes of ORCTs of *V. major*. The RMSEadj values calculated by Model 1 were smaller than those calculated by Model 2 for 346 out of 360 (ca. 96%) ORCT samples, and for the remaining 14 samples, the RMSEadj values between the two models were the same. The mean ln(RMSEadj) of Model 1, *M*1, was significantly smaller than that of Model 2, *M*<sup>1</sup> (*t* = −5.5666, df = 718, *p* < 0.001). This indicates that Model 1 with four parameters provided a better fit than Model 2 with three parameters, although the difference was significant in statistics but fairly small. The percentage error, i.e., |(*M*<sup>2</sup> − *M*1)/*M*1| × 100%, between the two mean ln(RMSEadj) values, was less than 2.7% (Figure 3b). From the tradeoff between the model complexity and goodness of fit, Model 1 is not as good as Model 2, and the latter has a more concise model structure.

**Figure 2.** Representative images of *V. major* flowers with five petals (top view; **a**,**b**), and the measured (gray) and predicted (red) outer rims of corolla tubes (ORCTs) using Model 1 (**c**,**d**) and Model 2 (**e**,**f**). The boundary coordinate data of ORCTs were obtained from scanned images, and were fitted by the R package 'biogeom'. RMSEadj is the adjusted root-mean-square error (Equation (7)).

**Figure 3.** Comparison of adjusted root-mean-square errors (**a**, RMSEadj, Equation (7)) and their ln-transformations between Models 1 and 2 (**b**). The thick horizontal lines represent median values in the boxes; a box's body length represents the difference between the 3/4 quantile and the 1/4 quantile; whiskers represent 1.5 times the box's body length or maximum (or minimum) values; and the small gray open circles represent the distribution of data points. The ln-transformed values in (**b**) were compared by a paired sample *t*-test.

The results of the linear regression of the data of ln(*n*ˆ <sup>3</sup>) vs. ln(*n*ˆ <sup>2</sup>) and the distribution diagrams of 2000 bootstrap replicates of the regression intercept (α) and slope (β) indicated that the 95% CI of the intercept included zero and the 95% CI of the slope included 1.0 (Figures 4 and 5). Thus, there was no significant difference between the estimated values of *n*<sup>3</sup> and that of *n*2, suggesting that the three-parameter Model 2 is sufficient to depict the geometry of the ORCT of *V. major*, although the goodness of fit of Model 1 is slightly greater than that of Model 2 (Figure 3). The estimated values of parameters in Models 1 and 2 have been listed in Tables S1 and S2 in the online Supplementary Materials.

**Figure 4.** Fitted results to the data of the estimated values of the model parameter *n*<sup>3</sup> vs. the model parameter *n*<sup>2</sup> (Equation (3)). The data were fitted according to reduced major axis protocols on a log-log scale. *y* represents the ln-transformation of the estimated value of *n*<sup>3</sup> ; *x* represents the ln-transformation of the estimated value of *n*<sup>2</sup> ; the straight line is the regression line; CIintercept represents the 95% confidence interval of the intercept; CIslope represents the 95% confidence interval of the slope; *r* 2 is the coefficient of determination that reflects the goodness of fit; and *N* is the sample size, i.e., the number of flowers used.

**Figure 5.** Frequency histograms of bootstrap replicates of the intercept (**a**) and slope (**b**) of the regression between the model parameters *n*<sup>3</sup> and *n*<sup>2</sup> (Equation (3); Figure 4). The red curves are normal density curves; and the blue vertical straight lines represent 0.025 quantile (left) and 0.975 quantile (right) that are the lower and upper bounds of the 95% confidence interval in each panel.

Although all *p* values of the normality test were smaller than 0.05, indicating that the distributions of residuals between the observed and predicted radii were not normal, the skewness only ranged from −1.1 to 1.1, which suggested that the distributions of residuals were not seriously skewed. The means of residuals for all samples ranged from −0.0002 to 0.0006 approximate to zero, and all of the 95% CIs of residuals included 0 (see Table S2). In other words, the distribution of residuals is approximate to be normal, and the above results at least demonstrated the validity of nonlinear least-squares in fitting the ORCTs.

### **4. Discussion**

### *4.1. Analysis of the Prediction Errors*

Multiple artificial measurement errors occurred when we scanned and extracted the boundaries of the 360 outer rims of corolla tubes (ORCTs). On the way from the site of plant growth to the laboratory, artificial physical pressures might have caused a certain deformation of the flowers. In addition, during scanning, the physical pressures on different contact points of a ORCT with the scanning plane could also cause a certain morphological deformation of the ORCT. Thus, there are some parts on an ORCT especially close to the vertices of the pentagon that are not well resolved in scanned images. This generates some uncertainties in extracting the ORCT. Nevertheless, given that such deformations were smaller than 0.05 mm, we conclude that the impact of such deformations on model fitting can be neglected.

In addition, the interspecific variation from a perfect symmetrical pattern in the growth process of *V. major* probably increased the deviation of the actual shape of ORCT from the predicted shape. The intraspecific variation reflects both genotypic differences as well as plastic modifications to differences in the growth microenvironment, including differences in the availability in light, water and nutrients [30–34]. Although the plants used were irrigated by sprinklers, the individuals closer to the water source formed a denser cover, indicating that water availability introduced some changes in the plant phenotype in this study. A stronger intraspecific competition for space at greater water availability caused

some physical compression within adjacent individuals during flower growth; this resulted in the deformation of some flowers to a certain degree and subsequently led ORCTs to be less symmetrical. All of the factors outlined likely had some impact on the prediction errors using the two simplified Gielis models to describe the shapes of ORCTs. Nevertheless, the RMSEadj values of the two models were all smaller than 0.05. This verified the validity of the Gielis equation in fitting the real natural geometries of corolla tubes of *V. major*.

### *4.2. Comparison of the Two Models*

The parameters *n*<sup>2</sup> and *n*<sup>3</sup> determine the extent of symmetry for the geometry generated by Model 1. The regression analysis indicated that there was no significant difference between *n*<sup>2</sup> and *n*3, indicating that the analysis can be simplified by replacing *n*<sup>3</sup> with *n*2. Therefore, we conclude that Model 2 with fewer parameters can replace Model 1 for describing the ORCT of *V. major*. Although most of the adjusted RMSE values calculated by Model 1 were slightly smaller than those calculated by Model 2 among 360 ORCTs, all of the RMSEadj values for both models were less than 0.05. This means that for any sample, the average absolute deviation between the observed and predicted radii from the polar point to an ORCT was less than 5% of the radius of a circle whose hypothetical area equals that of the ORCT.

The use of Model 2 has several practical advantages. In particular, as fewer parameters need to be fitted, the efficiency of parameter estimation is higher. This implies that the probability for the model parameters not to converge is lower for Model 2; also, the average computer running time to complete the optimization for each dataset using Model 2 is less than that using Model 1. In addition, Model 1 with an additional parameter *n*<sup>3</sup> tends to lead to an overfitting of the planar coordinate data of ORCTs. Accordingly, this blurs the biological meaning of the parameter values. Model 2 with parameters *n*<sup>3</sup> = *n*<sup>2</sup> diminishes the possibility of generating an asymmetrical curve and therefore agrees better with previous studies (e.g., refs. [7,9]).

### *4.3. Variation in the Number of Polygon Sides of the ORCT of V. major*

In addition to pentagonal ORCTs, there was a very small proportion of quadrangular ORCTs among the flowers of the sampled *V. major* population (Figure 6a). In addition, flowers with hexagonal ORCTs can be observed in this species (personal observations), but were not found in the current study. The plants with pentagonal ORCTs have five petals, while the plants with quadrilateral ORCTs have four petals; each angle of an ORCT corresponds to a petal (Figure 1and Figure 6). In flower sampling, we found only six specimens with quadrilateral petals, i.e., a ratio of 6:366 in the total sample of flowers collected. It is currently unclear what causes the formation of four-petal flowers in *V. major*, but given that it is a clonal species, it is unlikely that the variation of petal number is of genetic origin, although we cannot rule out that the sampled population was a mix of several clones. In fact, the situation is similar to *Syringa* spp. that typically have four-petal flowers, but occasionally form five-petal flowers, and sometimes even flowers with more than five petals in the same inflorescence.

Similar to flowers with pentagonal ORCTs, we fitted the coordinate data of the six quadrilateral ORCTs of *V. major* using Model 2. As with flowers with pentagonal ORCTs, all of the adjusted RMSE values for flowers with quadrilateral ORCTs were smaller than 0.05 (Table S3 in the online Supplementary Materials). Thus, Model 2 also well describes the shape of the ORCT of *V. major* with four petals (Figure 6b).

**Figure 6.** The vertical projection of *V. major* flower with four petals (**a**), and the comparison between the observed and predicted ORCTs (**b**). In panel (**b**), the gray curve represents the observed ORCT and the red curve represents the predicted ORCT using Model 2. Presentation and symbols as in Figure 3.

### **5. Conclusions**

The simplified Gielis equations with *m* = 5 or *m* = 4 are both valid in describing the outer rims of the corolla tubes (in top view) (ORCTs) of *V. major*. A further simplified Gielis model with *n*<sup>2</sup> = *n*<sup>3</sup> predicted a perfectly axial symmetrical ORCT. The use of the simpler model can avoid model overfitting and reduce the running time for completing parameter estimation. This work shows that the ORCTs of plants can be modeled by the Gielis equation and can provide a reference for future research on superelliptic shapes of plant organs or tissues, e.g., projections of calyxes and polygonal cross-sections of fruits, e.g., ref. [9]). Additionally, the small difference between the parameters *n*<sup>2</sup> and *n*<sup>3</sup> in Model 1 did not cause a large deviation for ORCTs from a perfectly axial symmetrical pattern. This work provides evidence that the natural morphology of plant organs and tissues follows the superformula proposed by Gielis [2], and the process of growth and development of corolla tubes exhibits a general but unknown biophysical mechanism that is commonly shared by other organs or tissues of many species, which at least can be reflected by the same mathematical model. Thus, it is worth further examining the validity of the Gielis equation in depicting more geometries found in nature to uncover the influence(s) of such a general biophysical mechanism on biological geometry and morphology.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants11151987/s1; Table S1: Fitted results of the planar coordinate data of ORCTs of *V. major* with five petals using Model 1; Table S2: Fitted results of the planar coordinate data of ORCTs of *V. major* with five petals using Model 2; and Table S3: Fitted results of the planar coordinate data of ORCTs of *V. major* with four petals using Model 2.

**Author Contributions:** Investigation, L.W. and Q.M.; formal analysis, L.W., Ü.N. and P.S.; writing—original draft preparation, L.W. and P.S.; writing—review and editing, Ü.N. and J.G.; supervision, Ü.N., J.G. and P.S.; L.W. and Q.M. contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data can be found in Tables S1–S3 in the online Supplementary Materials.

**Acknowledgments:** The authors would like to thank Yabing Jiao and Long Chen for their valuable help in flower sampling and data acquisition. The authors are also thankful to three reviewers for their extremely constructive comments for improving the quality of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


## *Article* **The Modified Brière Equation and Its Applications**

**Jun Jin <sup>1</sup> , Brady K. Quinn <sup>2</sup> and Peijian Shi 3,\***


**Abstract:** The Brière equation (BE) is widely used to describe the effect of temperature on the development rate of insects, and it can produce both symmetrical and asymmetrical bell-shaped curves. Because of its elasticity in curve fitting, the integrated form of BE has been recommended for use as a sigmoid growth equation to describe the increase in plant biomass with time. However, the start time of growth predicted by the sigmoid growth equation based on the BE is not completely comparable to empirical crop growth data. In the present study, we modified the BE by adding an additional parameter to further increase its elasticity for data fitting. We termed this new equation the modified Brière equation (MBE). Data for the actual height and biomass of 15 species of plants (with two cultivars for one species) were fit with the sigmoid growth equations based on both the BE and MBE assuming that the growth start time was zero for both. The goodness of fit of the BE and MBE sigmoid growth equations were compared based on their root-mean-square errors and the corresponding absolute percentage error between them when fit to these data. For most species, we found that the MBE sigmoid growth equation achieved a better goodness of fit than the BE sigmoid growth equation. This work provides a useful tool for quantifying the ontogenetic or population growth of plants.

**Keywords:** axial symmetry; curve fitting; ontogenetic growth; sigmoid curve; symmetry

### **1. Introduction**

The ontogenetic growth trajectories of animals and plants usually exhibit a sigmoid pattern, and many mathematical equations have been proposed to describe the changes in growth data (e.g., biomass, diameter at breast height, height, length, etc.) with time [1–6]. Among these, the logistic equation is perhaps the most commonly used [5–8]. The threeparameter logistic equation assumes that the growth rate versus time curve is perfectly symmetrical but is not suited to all empirical datasets [9]. Thus, some equations predicting asymmetrical (skewed) growth rate curves have been used instead to reflect the growth trajectories of animals and plants [3,9,10].

In thermal biology, the development (or growth) rate (i.e., the proportion or amount of development or growth, respectively, completed per unit time) of organisms has been demonstrated to be an asymmetrical (usually left-skewed), bell-shaped function of temperature [11,12]. There are many different mathematical models that can produce skewed and symmetrical development (or growth) rate versus temperature curves [13,14]. Among these temperature-dependent development rate models, the equation proposed by Brière et al. [15] is relatively simple, yet it can produce an asymmetrical curve, includes biologically meaningful parameters, and is able to describe the effect of temperature on the development rate of many organisms:

$$r(T) = a\mathbf{x}(T - T\_{\text{min}})(T\_{\text{max}} - T)^{1/m},\tag{1}$$

**Citation:** Jin, J.; Quinn, B.K.; Shi, P. The Modified Brière Equation and Its Applications. *Plants* **2022**, *11*, 1769. https://doi.org/10.3390/plants11131769

Academic Editor: Oleg Chertov

Received: 15 June 2022 Accepted: 1 July 2022 Published: 3 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where *r*(*T*) represents the development rate at temperature *T*; *T*min and *T*max represent the lower and upper threshold temperatures of development, respectively, and *r*(*T*) = 0 when *T* < *T*min or *T* > *T*max; and *a* and *m* are parameters to be estimated. Brière et al. [15] suggested that a simplified version of Equation (1) could be used by fixing *m* = 2; this simplified equation, referred to as the Brière-1 equation (versus the original Brière equation, or Brière-2 equation), is not considered further herein. All subsequent references to the Brière equation hereinafter refer to the original Brière equation, or Brière-2 equation. Up to 25 May 2022, the original paper by Brière et al. [15] has been cited 745 times according to Google Scholar.

Yin et al. [3] introduced a new sigmoid growth equation using the integrated form of the beta equation while replacing temperature with time as the independent variable, and it has been shown to provide a good fit to the actual growth data of plants [3,9,16,17]. Similarly, Cao et al. [18] presented a sigmoid growth equation based on the Brière equation. However, although the new growth equation introduced by Cao et al. [18] could fit the growth data of some crops very well, the start times of growth predicted by it were not reliable as the predicted start time seriously deviated from the actual sowing time [18]. Further, when the start time is manually set to zero, the prediction errors for the sigmoid equation proposed by Cao et al. [18] relative to empirical data can be quite large.

In this study, we propose a modification of the Brière equation by adding an extra parameter to render the new equation more flexible in curve fitting and reduce the dependency of the curve fitting on the start time of growth. We used empirical data from 15 species of plants (with two cultivars for one species) to test the validity of a new sigmoid growth equation based on this modified Brière equation (MBE) and compared the goodness of fit of the sigmoid growth equation based on the MBE to that based on the original Brière equation (BE).

### **2. Materials and Methods**

### *2.1. Plant Materials*

The seeds of 11 species of common agricultural crops (with two cultivars for one species, 12 crops in total), including sunflower (*Helianthus annuus*), peanut (*Arachis hypogaea*), black soybean (*Glycine max* 'Kuromame'), soybean (*Glycine max*), kidney bean (*Phaseolus vulgaris*), garden pea (*Pisum sativum*), adzuki bean (*Vigna angularis*), mung bean (*Vigna radiata*), cotton (*Gossypium herbaceum*), sweet sorghum (*Sorghum bicolor*), corn (*Zea mays*), and Mexican corn (*Zea mexicana*), in northern China were sown in a field at Jinan, China on 27 June 2011. We then randomly sampled 20 individuals from each species (or each cultivar) to measure their dry mass on each of 15 subsequent investigation dates from 11 July to 20 September 2011 once every five days on average. In total, 3000 individual samples were collected. All individuals were sampled between 7:00 a.m. to 8:00 a.m., and the roots were washed with running water to remove soil. For small crops, whole plants were dried for 24–48 h at 60 ◦C. However, for large crops, such as sunflower, sweet sorghum, corn, and Mexican corn, whole plants were dried for 72 h at 80 ◦C. Here, we defined the crops with mean whole-plant dry mass ≥ 150 g and mean aboveground height ≥ 1.5 m on the last investigation date (i.e., 85 days from the sowing time) as large crops and those with mean whole-plant dry mass < 150 g and mean aboveground height < 1.5 m at this time as small crops. Samples were weighed after drying to obtain their dry mass. Shi et al. [5] can be consulted for further details of the sampling procedures.

We also used the height data from four species of bamboo (*Phyllostachys iridescens*, *Phyllostachys mannii*, *Pleioblastus maculatus*, and *Sinobambusa tootsik*) previously collected at different investigation dates [9]. The bamboo shoots were grown at the Nanjing Forestry University campus in the spring of 2016. We defined the time when a bamboo shoot tip was first observed at ground level as time = 1, and the previous day as time = 0. We measured the height of each shoot at 12:00 p.m. every day at the early growth stage, and then measured the height at 12:00 p.m. once every three days when the height changed more slowly.

The above datasets were obtained from the 'crops' and 'shoots' datasets, respectively, included in the R (version 4.2.0) [19] package 'IPEC' (version 1.0.3) [20]. The above datasets were obtained from the 'crops' and 'shoots' datasets,respectively, included in the R (version 4.2.0) [19] package 'IPEC' (version 1.0.3) [20].

was first observed at ground level as time = 1, and the previous day as time = 0. We meas‐ ured the height of each shoot at 12:00 p.m. every day at the early growth stage, and then measured the height at 12:00 p.m. once every three days when the height changed more

*Plants* **2022**, *11*, x FOR PEER REVIEW 3 of 10

#### *2.2. Methods 2.2. Methods*

slowly.

The original Brière equation was proposed to describe the effect of temperature on the development rate of many organisms, especially arthropods [13,15,21] (see Equation (1)). After replacing temperature *T* with time *x*, and the development rate *r*(*T*) with the growth rate *f*(*x*) in the original Brière equation, we have: The original Brière equation was proposed to describe the effect of temperature on the development rate of many organisms, especially arthropods [13,15,21] (see Equation (1)). After replacing temperature *T* with time *x*, and the development rate *r*(*T*) with the growth rate *f*(*x*) in the original Brière equation, we have:

$$f(\mathbf{x}) = a\mathbf{x}(\mathbf{x} - \mathbf{x}\_{\text{min}})(\mathbf{x}\_{\text{max}} - \mathbf{x})^{1/m},\tag{2}$$

where *f*(*x*) represents the growth rate at time *x*; *x*min and *x*max represent the start and end times of growth, respectively, and *f*(*x*) = 0 when *x* < *x*min or *x* > *x*max; and *a* and *m* are parameters to be estimated. where *f*(*x*) represents the growth rate at time *x*; *x*min and *x*max represent the start and end times of growth, respectively, and *f*(*x*) = 0 when *x* < *x*min or *x* > *x*max; and *a* and *m* are param‐ eters to be estimated.

Herein, we propose including an additional parameter to be estimated, δ, in Equation (2) to render it more elastic in curve fitting as follows: Herein, we propose including an additional parameterto be estimated, δ, in Equation (2) to render it more elastic in curve fitting as follows:

$$f(\mathbf{x}) = a \left| \mathbf{x} (\mathbf{x} - \mathbf{x}\_{\min}) (\mathbf{x}\_{\max} - \mathbf{x})^{1/m} \right|^{\delta}. \tag{3}$$

We refer to Equations (2) and (3) as the BE and MBE for convenience hereinafter. The integrated forms of the BE and MBE were used to fit the dry mass versus time data of crops and those of height versus time data of bamboo shoots. Figure 1A illustrates the influence of the new δ parameter on the curve shapes plotted by the modified BE. We refer to Equations (2) and (3) as the BE and MBE for convenience hereinafter. The integrated forms of the BE and MBE were used to fit the dry mass versus time data of crops and those of height versus time data of bamboo shoots. Figure 1A illustrates the influence of the new δ parameter on the curve shapes plotted by the modified BE.

**Figure 1.** Influence of the new additional parameter, δ, on the shape of the curves plotting the mod‐ ified Brière equation (**A**), and on those plotting the sigmoid equation based on the integrated form of the modified Brière equation (**B**). Here, *a* = 0.01, *m* = 3, *x*min = 0, and *x*max = 35. In this example, no specific units are used; the *x*‐axis represents a generic time variable (hours, days, weeks, etc.) and the *y*‐axes generic (**A**) growth rate (cm day<sup>−</sup>1, g h<sup>−</sup>1, etc.) and growth (cm, g, etc.) metrics. **Figure 1.** Influence of the new additional parameter, δ, on the shape of the curves plotting the modified Brière equation (**A**), and on those plotting the sigmoid equation based on the integrated form of the modified Brière equation (**B**). Here, *a* = 0.01, *m* = 3, *x*min = 0, and *x*max= 35. In this example, no specific units are used; the *x*-axis represents a generic time variable (hours, days, weeks, etc.) and the *y*-axes generic (**A**) growth rate (cm day−<sup>1</sup> , g h−<sup>1</sup> , etc.) and growth (cm, g, etc.) metrics.

After integrating the BE and MBE, two sigmoid growth equations are obtained, i.e., After integrating the BE and MBE, two sigmoid growth equations are obtained, i.e.,

$$y = \begin{cases} 0 & \text{x < x\_{\text{min}}}\\ \int\_{\text{x\_{\text{min}}}}^{\text{x}} f(\mathbf{x}) \, d\mathbf{x} & \text{x \in [x\_{\text{min}}, x\_{\text{max}}]} \,. \end{cases} \tag{4}$$

*x x*

 max Here, we obtained two sigmoid growth equations when using Equations (2) and (3) to represent *f*(*x*), respectively. There is an explicit mathematical expression for the sigmoid equation based on the integrated form of BE [3,17], which is referred to as the BE sigmoid equation hereinafter for convenience. However, there is no explicit mathematical form for the integrated form of MBE. Thus, we calculated the numerical integral of MBE, which we refer to as the MBE sigmoid equation for convenience. Figure 1B illustrates the curves produced by the MBE sigmoid equation with different δ values.

The Nelder–Mead optimization method [22] was used to fit the above equations to each plant dataset by minimizing the residual sum of squares (RSS) between observed and predicted *y* values. For the empirical crop biomass and bamboo shoot height data we used, the start time of growth is actually known. Thus, we fixed the starting time to be zero, in which case Equation (3) had four parameters and Equation (2) had three parameters. We then calculated the root-mean-square error (RMSE) for each dataset:

$$\text{RMSE} = \sqrt{\text{RSS}/n} \tag{5}$$

where *n* represents the number of data points. To examine whether the additional parameter in Equation (3) was warranted relative to using the simpler Equation (2), we used the percentage error of the absolute difference (APE, in %) between the RMSE of the BE sigmoid equation and that of the MBE sigmoid equation [23,24] as follows:

$$\text{APE} = \frac{|\text{RMSE}\_{\text{BE}} - \text{RMSE}\_{\text{MBE}}|}{\text{RMSE}\_{\text{BE}}} \times 100\%. \tag{6}$$

As a rule of thumb, when APE > 5%, the introduction of the additional parameter (δ) in Equation (3) versus (2) is justified; otherwise, it is not worthwhile.

The whole-plant dry mass tends to be proportional to the whole-plant fresh mass [5,25] if there is no water loss. In that case, if the sigmoid equation based on the integrated form of MBE can fit the dry mass data well, it should also be applicable to fresh mass data. The difference here would only affect the numerical value of the parameter *a* in Equation (3). However, given the water loss during crop sampling, the dry mass is more reliable than the fresh mass. Thus, for the crops, we used the dry mass rather than the fresh mass to test the validities of the BE and MBE sigmoid equations.

The package 'biogeom' (version 1.0.5) [26,27] was used in R (version 4.2.0) [19] to estimate model parameters, and R (version 4.2.0) was also used to carry out all other calculations.

### **3. Results**

For each of the crop species tested except kidney bean (11/12 datasets), the MBE sigmoid equation had a lower RMSE than the BE sigmoid equation, and the APE was greater than 5% (Figure 2). This confirmed the validity of the MBE sigmoid equation and its superiority relative to the BE sigmoid equation. Because of the lack of biomass data at the mature stages of sunflower and peanut, both sigmoid functions overestimated the end time of growth, leading the predicted curves to exhibit exponential growth at the early and middle growth stages (Figure 2A,B).

For each of the four species of bamboo considered, the MBE sigmoid equation also had a lower RMSE than the BE sigmoid equation, and the APE was larger than 5% for three of the four species (Figure 3). For *Pleioblastus maculatus*, the APE was equal to 4.68%, or approximately 5%. These results demonstrate that the BE and MBE sigmoid equations are both suited to representing data for growth in terms of height for bamboo shoots, at least provided that the initial time at which growth starts (i.e., *x* = 0) can be accurately known. Nevertheless, the MBE sigmoid equation is still better.

The estimated values of the parameters of the BE sigmoid equation and the MBE sigmoid equation for the 15 species of plants (with two cultivars for one species) tested are listed in Tables S1 and S2 in the online Supplementary Materials.

**Figure 2.** Results of fitting two sigmoid functions (based on the BE and MBE) to the whole‐plant dry mass versus time data of 11 crop species (with two cultivars for one crop species). The small open circles represent the observed dry mass at different times after the sowing date; the red dashed curves represent the dry mass values predicted by the BE sigmoid equation; the blue solid curves represent the dry mass values predicted by the MBE sigmoid equation. RMSE represents the root‐ mean‐square error between the observed and predicted *y* values; APE represents the percentage **Figure 2.** Results of fitting two sigmoid functions (based on the BE and MBE) to the whole-plant dry mass versus time data of 11 crop species (with two cultivars for one crop species). The small open circles represent the observed dry mass at different times after the sowing date; the red dashed curves represent the dry mass values predicted by the BE sigmoid equation; the blue solid curves represent the dry mass values predicted by the MBE sigmoid equation. RMSE represents the root-mean-square error between the observed and predicted *y* values; APE represents the percentage error of the absolute difference between the two equations' RMSE values; *n* represents the sample size. Panels (**A**–**L**) represent different crops.

error of the absolute difference between the two equations' RMSE values; *n* represents the sample

size. Panels (**A**–**L**) represent different crops.

**Figure 3.** Results of fitting two sigmoid functions (based on the BE and MBE) to the shoot height versus time data of four species of bamboo. The small open circles represent the observed shoot height at different times (days) after the date when the shoot tip first emerged from the soil; the red dashed curves represent the height values predicted by the BE sigmoid equation; the blue solid curves represent the height values predicted by the MBE sigmoid equation. RMSE represents the root‐mean‐square error between the observed and predicted *y* values; APE represents the absolute percent error difference between the two equations' RMSE values; *n* represents the sample size. Panels (**A**–**D**) represent different bamboo species. **Figure 3.** Results of fitting two sigmoid functions (based on the BE and MBE) to the shoot height versus time data of four species of bamboo. The small open circles represent the observed shoot height at different times (days) after the date when the shoot tip first emerged from the soil; the red dashed curves represent the height values predicted by the BE sigmoid equation; the blue solid curves represent the height values predicted by the MBE sigmoid equation. RMSE represents the root-mean-square error between the observed and predicted *y* values; APE represents the absolute percent error difference between the two equations' RMSE values; *n* represents the sample size. Panels (**A**–**D**) represent different bamboo species.

#### **4. Discussion 4. Discussion**

#### *4.1. Elasticity in Curve Fitting of the Two Sigmoid Equations 4.1. Elasticity in Curve Fitting of the Two Sigmoid Equations*

The inclusion of an additional parameter, δ, in the MBE greatly increased the elastic‐ ity of the MBE sigmoid equation in curve fitting; i.e., it increased the range of sigmoid curves with different curvatures that can be fit relative to the BE sigmoid equation (Figure 1). The larger the numerical value of δ is, the greater the curvature of the line will be that is generated by the MBE sigmoid equation. For 14 of the 16 tested datasets, the APE values were greaterthan 5%, which indicates that the increased elasticity in curve fitting achieved by the additional parameter in the MBE sigmoid equation outweighed the cost of the model's increased complexity [28]. In fact, the elasticities of the two sigmoid equations in curve fitting mainly depend on the abilities of the BE and MBE to accurately describe growth rates. However, it is more difficult to directly measure the growth rate per unit time in practice than to measure accumulated biomass or height after some time interval has elapsed. Although Cao et al. [18] showed that the BE sigmoid growth function achieved a good fit to the dry mass data of six of the twelve datasets of crops used herein, The inclusion of an additional parameter, δ, in the MBE greatly increased the elasticity of the MBE sigmoid equation in curve fitting; i.e., it increased the range of sigmoid curves with different curvatures that can be fit relative to the BE sigmoid equation (Figure 1). The larger the numerical value of δ is, the greater the curvature of the line will be that is generated by the MBE sigmoid equation. For 14 of the 16 tested datasets, the APE values were greater than 5%, which indicates that the increased elasticity in curve fitting achieved by the additional parameter in the MBE sigmoid equation outweighed the cost of the model's increased complexity [28]. In fact, the elasticities of the two sigmoid equations in curve fitting mainly depend on the abilities of the BE and MBE to accurately describe growth rates. However, it is more difficult to directly measure the growth rate per unit time in practice than to measure accumulated biomass or height after some time interval has elapsed. Although Cao et al. [18] showed that the BE sigmoid growth function achieved a good fit to the dry mass data of six of the twelve datasets of crops used herein, the estimated growth start times were later than the actual sowing time and the times seedings appeared (see Figure 1 in Ref. [18]). It is apparently reasonable to fix the start time to the sowing time given that all the seeds were planted on the same day. In our study, the MBE sigmoid

function achieved a good fit to observations without estimating unreasonable start times, and it showed good elasticity in curve fitting through changes in the numerical value of the δ parameter. After all, the *x*min parameter of the BE and MBE has an explicit geometrical and biological meaning, so it is unreasonable to over- or underestimate its numerical value. An over- or underestimated *x*min value can easily mislead the user of such equations to predict inaccurate growth. On the other hand, the new δ parameter has no clear biological or geometrical meaning, and it only serves as a constant to be estimated to improve the flexibility of the equation in curve fitting.

### *4.2. Why We Did Not Compare the Two Sigmoid Equations with Other Equations*

As stated in the Introduction section, there are many temperature-dependent development (or growth rate) equations that can produce a skewed curve [9,13,14,21]. Shi et al. [17] found that the sigmoid equation based on the integral of the beta equation [3] can achieve better goodness of fit than traditional growth equations, including the exponential and logistic equations [8], Gompertz equation [29], von Bertalanffy equation [1], and ontogenetic growth equation [2]. Shi et al. [9] demonstrated that the sigmoid equation based on the integral of the modified beta equation and that of the modified Lobry–Rosso–Flandrois (LRF) equation [30,31] can be better-suited to representing the ontogenetic growth data of animals and plants. However, the mathematical expressions of the BE and MBE are simpler than those of the beta equation, modified beta equation, LRF equation, and modified LRF equation (Equations (2) and (3) versus those published in Ref. [9]). In the present study, we did not compare the BE and MBE sigmoid equations with other sigmoid equations because we mainly wanted to examine the strengths and drawbacks of the BE and MBE sigmoid equations and to test whether the latter is better than the former for curve fitting. A systematic comparison among sigmoid equations is worth carrying out using more datasets in future studies.

### *4.3. Reliability of Estimated Parameters in the Two Sigmoid Equations*

Although the estimated parameters, especially the start and end times of growth, are considered to be meaningful for predicting when these time points will occur during the growth of plants, the reliability of the prediction is largely constrained by whether the dataset includes a full range of growth data, especially those at the mature stage (i.e., the asymptotic values of plant biomass or height). If the asymptotic values of plant size at the mature stage are known, the predicted end time of growth tends to be reliable [6,9,18]; however, if such data are lacking, the predicted end time tends to be overestimated (see Figure 2A,B,L). To reduce the uncertainty of predictions and the complexity of the model, we suggest fixing the starting time and allowing the δ parameter to vary. It is more important to estimate the end time rather than the start time of growth for a crop product or forest resource management. If the estimated end time of growth is reliable, it can produce a reliable growth rate curve with which one can predict when plants reach their maximum growth rate.

### *4.4. Other Potential Applications of the Modified Brière Equation*

Given that the MBE is flexible in fitting a skewed bell-shaped curve, it is potentially useful to model the profiles of ovate leaves [32]. Shi et al. [32] formed two axially symmetrical curves using the modified beta equation (also the modified LRF equation) to fit the ovate leaf shapes of *Neocinnamomum* plants (Lauraceae). Given that the MBE can generate similar skewed curves, the MBE should be similarly useful in describing ovate leaf shape. Relative to the BE, the MBE is more flexible due to its inclusion of an additional parameter, δ, as discussed above, and it can thus also fit the actual boundary coordinates of ovate leaves (Figure 4). Like the original BE, the MBE may also be useful in data-fitting and modelling of temperature-dependent development and growth rates in animals and plants [15,21]. While such additional applications of the MBE are beyond the scope of the present work (and, thus, we do not discuss them in detail here), it will be worthwhile for

future studies to assess the validity of the MBE in fitting ovate and obovate leaf shapes and temperature-dependent development rate data. future studies to assess the validity of the MBE in fitting ovate and obovate leaf shapes and temperature‐dependent development rate data.

leaves (Figure 4). Like the original BE, the MBE may also be useful in data‐fitting and modelling of temperature‐dependent development and growth rates in animals and plants [15,21]. While such additional applications of the MBE are beyond the scope of the present work (and, thus, we do not discuss them in detail here), it will be worthwhile for

*Plants* **2022**, *11*, x FOR PEER REVIEW 8 of 10

**Figure 4.** Results of fitting the MBE to the boundary coordinates of the leaves of four species of *Neocinnamomum*. The gray curves are scanned (actual) leaf perimeters, and the red curves are leaf perimeters predicted by the MBE. The boundary coordinate data came from the dataset '*Neocin‐ namomum*' in R package 'biogeom' (https://cran.r‐project.org/web/packages/biogeom/index.html; accessed on 30 May 2022). **Figure 4.** Results of fitting the MBE to the boundary coordinates of the leaves of four species of *Neocinnamomum*. The gray curves are scanned (actual) leaf perimeters, and the red curves are leaf perimeters predicted by the MBE. The boundary coordinate data came from the dataset '*Neocinnamomum*' in R package 'biogeom' (https://cran.r-project.org/web/packages/biogeom/index.html; accessed on 30 May 2022).

#### **5. Conclusions 5. Conclusions**

A better goodness of fit was obtained with the MBE sigmoid equation for 15 of the 16 datasets of plants than the BE sigmoid equation, and it had the same RMSE as that of the BE sigmoid equation for the one remaining species. For most datasets (14/16), the percent‐ age errors of the absolute difference between the RMSE of the MBE sigmoid equation and that of the BE sigmoid equation were greater than 5%, which indicates that the addition of the δ parameter to the original Brière equation was worthwhile, improving the validity of the MBE sigmoid equation in reflecting the actual growth data of plants. In addition, by virtue of the estimated model parameters, whether the growth curve of a given species is left‐ or right‐skewed can be evaluated (i.e., based on the derivative of the MBE sigmoid equation), and, thus, one can use this to predict the time associated with the maximum A better goodness of fit was obtained with the MBE sigmoid equation for 15 of the 16 datasets of plants than the BE sigmoid equation, and it had the same RMSE as that of the BE sigmoid equation for the one remaining species. For most datasets (14/16), the percentage errors of the absolute difference between the RMSE of the MBE sigmoid equation and that of the BE sigmoid equation were greater than 5%, which indicates that the addition of the δ parameter to the original Brière equation was worthwhile, improving the validity of the MBE sigmoid equation in reflecting the actual growth data of plants. In addition, by virtue of the estimated model parameters, whether the growth curve of a given species is left- or right-skewed can be evaluated (i.e., based on the derivative of the MBE sigmoid equation), and, thus, one can use this to predict the time associated with the maximum growth rate.

growth rate. **Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants11131769/s1, Table S1: The estimates of the model parameters for the whole-plant dry mass versus time data of 12 crops; Table S2: The estimates of the model parameters for the shoot height versus time data of four species of bamboo.

**Author Contributions:** Conceptualization, J.J., P.S. and B.K.Q.; methodology, P.S.; software, P.S. and B.K.Q.; formal analysis, J.J. and P.S.; investigation, J.J., P.S. and B.K.Q.; writing—original draft

preparation, J.J. and P.S.; writing—review and editing, B.K.Q.; funding acquisition, J.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key R&D Program of China (grant numbers 2020YFD1100405 and 2020YFD1100400).

**Data Availability Statement:** The data used in the present work have been packaged in package 'IPEC' (version 1.0.3; https://cran.r-project.org/web/packages/IPEC/index.html; accessed on 24 May 2022). The crop biomass growth data were tabulated in dataset 'crops', and the bamboo height growth data were tabulated in dataset 'shoots'.

**Acknowledgments:** We thank Meilin Fan, Shuiyuan Fang, Feng Ge, Xingyuan Men, Fang Ouyang, David A. Ratkowsky, and Chengkang Wang for their valuable help during the preparation of this work. We also thank the academic editor and four reviewers for their valuable comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland www.mdpi.com

*Plants* Editorial Office E-mail: plants@mdpi.com www.mdpi.com/journal/plants

Disclaimer/Publisher's Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Academic Open Access Publishing

mdpi.com ISBN 978-3-0365-9423-1