Indica rice (Oryza Sativa L.) cultivar “Yangdao 6” was used as the material. The rice seeds were germinated and then grew in a growth chamber under 16-h light, 28°C, and 70% humidity using hydroponics. The composition of the nutrient solution is determined following the method proposed by Yoshida et al. (1976). The pH value was adjusted to 5.8, and the dissolved oxygen concentration was maintained at 6.5–7.0 mg L^–1. Rice roots of 2-cm length from 20-day-old seedlings were marked using an oil-based marker pen. The marked roots were sampled for the anatomical observation after they grew to different lengths from 3 to 16 cm in the growth chamber.

Roots with different lengths were fixed in glutaraldehyde for 48 h and dehydrated in ethanol at the concentrations of 30, 40, 50, 60, and 70% for 20 min (each concentration). The cross-sections were prepared at every 0.5-cm interval from the root apex by hand sectioning with a razor blade under a pose microscope (Leica EZ24W, Germany) (Yamauchi et al., 2016). Each section was photographed using a microscope equipped with an AxioCam HRc camera (ZEISS, United States). All the images were analyzed by dividing the airspace area, and the total area of root sections and their areas were quantified with ImageJ (version 1.43) (Chen et al., 2021).

Four non-linear models were evaluated to describe the longitudinal pattern of aerenchyma formation. The mean values of the percentage of the aerenchyma area to the cross-sectional area (i.e., aerenchyma percentage) were fitted using cftool in MATLAB 2018a (MathWorks, Natick, MA, United States) with the following models:

where y is the dependent variable (i.e., aerenchyma percentage) at the distances from the root apex, x is the distance from the root apex, A is the parameter that represents the asymptotic value, B is the parameter related to the growth rate, and C is the parameter that represents the distance between the initial value and the asymptote (Tjørve and Tjørve, 2017; Sari et al., 2019; Teixeira et al., 2021).

With cftool in the MATLAB software, R² and the sum of squares of residuals (SSR) were obtained. The Akaike information criterion (AIC) and the Bayesian information criterion (BIC) were calculated as follows:

where R² is the coefficient of determination of the model, SSR is the sum of squares of residuals of the model, p is the number of the model parameters, AIC is the Akaike information criterion, and BIC is the Bayesian information criterion (Pinho et al., 2014).

The parameters of different models were estimated by the ordinary least squares (OLS) method and the bootstrap-based parameter estimation (Bates and Watts, 2007). In this procedure, the distribution plots of the obtained 1,000 bootstrap-based parameter estimates were constructed for each model, and the 95% CIs were computed as the difference between the 97.5th and 2.5th percentiles of 1,000 parameter estimates. The probability density function (short for density) is defined by the following formula: P(a≤x≤b)=∫abf(x)dx, where (a, b) is the interval in which x lies, P (a ≤ x ≤ b) is the probability that some value of x lies within this interval, dx is b−a, and f(x) is density.

The bias of the parameters was investigated by using the following: (1) the degree of symmetry of the sampling distribution of parameter estimates and (2) the difference between the average mean of the bootstrap parameter estimates and the OLS estimates. Asymmetric distributions and high differences between the mean bootstrap estimates and the OLS estimates were related to biased parameters (Ritz and Streibig, 2008; Sari et al., 2019).

The abscissa values (i.e., distance from the root apex) of the critical points in the Ti-Gompertz model were calculated by the following equations (Bem et al., 2017; Tjørve and Tjørve, 2017).

The initial points (InP) were obtained by equaling the dependent variance (i.e., aerenchyma percentage) to 1 as follows:

The maximum value points (MVP) were obtained by setting the dependent variance (i.e., aerenchyma percentage) to the 95% estimated parameter A as follows:

The inflection points (IPs) were set equal to the value of the parameter C.

The maximum acceleration points (MAPs) and the maximum deceleration points (MDPs) were obtained by equating the third-order derivative to zero as follows:

The experimental data were sorted using Microsoft Excel 2013. The models were fit and plotted using cftool in MATLAB software. The non-linear regression analysis such as the OLS and the bootstrap-based parameter estimations was conducted using the nlstools package in R 4.0.2, and the graphs were constructed using the ggplot2 package in R (Baty et al., 2015).

Root cross-sections were prepared every 0.5 cm from the apex of 2- to 16-cm-long adventitious roots (Figure 1A). According to all values of the percentages of the aerenchyma area to the cross-sectional area (short for aerenchyma percentage), the minimum and maximum values were 1.0 and 45.2%, respectively, as observed in Figures 1B,C. No aerenchyma was observed in roots ≤3 cm. The aerenchyma began to form in roots ≥4 cm. The InP of aerenchyma formation were 1.5 and 1.0 cm from the root apex of 4- to 10-cm roots and 11- to 16-cm roots, respectively. The aerenchyma percentage gradually increased as the distance from the root apex in roots that were 4–7 cm in length increased. For roots ≥8 cm, the aerenchyma percentage increased and then reached a maximum level (Figures 1A–C and Supplementary Figure 1). Comparing the aerenchyma percentage at the same point from the apex of roots with different lengths, the longer the root length was, the more developed the aerenchyma (Figure 1A). Additionally, the aerenchyma percentage was lower at approximately 1 cm from the basal part of all roots (Figure 1A).

Generally, the curves describing aerenchyma formation were S-shaped, except for the last 1 cm from the root apex (Supplementary Figure 1). The aerenchyma percentage data of the roots with different lengths were fit using four types of non-linear models such as two types of Gompertz (i.e., W₀-Gompertz and Ti-Gompertz), logistic, and von Bertalanffy models (Figure 2). The W₀-Gompertz model and the Ti-Gompertz model are the different reparameterization forms of the traditional Gompertz model. Our results showed that the fitted curves of the Ti-Gompertz and W₀-Gompertz models were the same and coincided with the observed values (Figure 2 and Supplementary Figure 1). A logistic model fits the aerenchyma data of roots <7 cm well (Figures 2A–C), but the fitted curves showed differences from the observed values at 0–1 cm from the apex of roots ≥7 cm. The maximum fitted aerenchyma percentage was 2.0%, while the observed value was 0 (i.e., no aerenchyma formation was observed) (Figures 2D–M). The von Bertalanffy model fits the aerenchyma data of roots ≤13 cm well (Figures 2A–J) but showed the maximum difference between fitted and observed values by 2.6% at 1–6 cm from the apex of roots >13 cm (Figures 2K–M).

Aiming to identify the best fit model, we used several procedures to check the goodness of fit and estimate the parameters. There was no difference between the W₀-Gompertz and Ti-Gompertz models based on the observed criteria. The R² values were greater than 0.99 (close to 1) for all models. Compared with logistic and von Bertalanffy models, two Gompertz models showed the lowest values for the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) within the same root length. Overall, there was a good adjustment of the investigated models, and the Gompertz models produced the best fit of the data for the aerenchyma percentage based on the three criteria (Supplementary Table 1). Furthermore, the difference in the estimates between the OLS and bootstrap-based approaches in the four models showed less bias for roots ≥6 cm (Figure 3). The bias of the estimates and the 95% CIs of the parameters for roots <6 cm among the four models were different (Figure 3). The W₀-Gompertz and von Bertalanffy models had larger biases and CIs than the other models (Figure 3). The Ti-Gompertz and logistic models had parameter estimates close to being unbiased, and all assumptions were met (Figure 4). Most of the distributions of the parameters were symmetrical except that parameter A in the von Bertalanffy model and parameter C in the W₀-Gompertz model were asymmetrical (Figure 4).

Consequently, the Ti-Gompertz model had higher goodness of fit, lower difference between parameter estimates by the OLS and bootstrapping, and symmetrical distribution of all parameters. This was a direct result of the good linear approximation of the Ti-Gompertz model compared to other models, suggesting that the Ti-Gompertz model was the best to describe the data of the aerenchyma formation.

The linear regression analysis was performed to unify the Ti-Gompertz model for roots with different lengths. Parameter A represents the asymptotic maximum aerenchyma percentage. For roots <13 cm, the values of the asymptotes increased as the root length (L) increased, and the equation was A = 13.02 + 2.54 L (P < 0.001; Figure 5A). The maximum aerenchyma percentage was 44.8 for all roots ≥13 cm and showed no correlation with the root length (A = 44.8 − 0.021 L; P = 0.934; Figure 5A). Parameter B represents the maximum speed of aerenchyma formation. For roots <13 cm, B had no correlation with the root length, and the maximum speed of aerenchyma formation was 0.751 (B = 0.751 + 0.0022 L; P = 0.222; Figure 5B). For roots ≥13 cm, B had a significantly positive correlation with the root length (B = −0.748 + 0.119 L; P < 0.001; Figure 5B). Parameter C represents the distance from the root apex when the aerenchyma formation speed reaches a maximum. C showed a significantly negative correlation with the root length, and the equations were C = 3.2 − 0.047 L (L < 13; P < 0.001) and C = 4.43 − 0.139 L (L ≥ 13; P < 0.001) (Figure 5C).

According to the minimum value of 1.0% for the aerenchyma percentage, we calculated the abscissa value of the initial point of aerenchyma formation under a dependent variance (i.e., aerenchyma percentage) of 1.0%. The coordinates of MAPs, IPs, and MDPs of the Ti-Gompertz model were obtained using the second- and third-order partial derivatives. The abscissa values of the MAPs, IPs, and MDPs appeared at 1.38–1.76 cm, 2.26–2.96 cm, and 3.13–4.19 cm from the root apex, respectively, and decreased as the root length increased (Table 1). The development of aerenchyma increased exponentially and logarithmically from the MAP to the IP and from the IP to the MDP, respectively (Table 1). These critical points in the Ti-Gompertz model suggested that the longer the roots were, the more developed the aerenchyma at the same point from the apex of roots with different lengths.

The entire period of aerenchyma formation was further divided into five stages as follows: the lag stage (Stage I), the starting stage (Stage II), the rapid stage (Stage III), the temperate stage (Stage IV), and the plateau (Stage V) (Figure 6). Overall, for adventitious roots <11 cm, there was no aerenchyma formation at <1.0 cm from the root apex (Figure 1). Stage II was defined as from the InP to the MAP, at which the aerenchyma percentage increased by 2.8–6.7% per cm. Stage III started at 1.38–1.76 cm from the root apex and continued to grow for approximately 2.0 cm, at which the aerenchyma percentage increased by 5.3–15.7% per cm. Stage IV was defined as from the MDP to the 95% MVP, at which the aerenchyma increased 3.3–6.8% per cm. For adventitious roots ≥8 cm, the aerenchyma formation reached a plateau at which point the aerenchyma percentage reached its maximum (Figure 2, Table 1, and Supplementary Figure 2). The initial point of the plateau was close to the root apex as the root length increased (Table 1). Furthermore, using the integral function of the Ti-Gompertz model, the volume of the aerenchyma was computed with the length and average diameter of the roots implying that once the root length and diameter were determined, the volume of the aerenchyma in the entire roots could be further computed (Figure 7).

Models are essential in biological science to summarize the data, to allow rapid predictions, and to design processes, such as “what-if” scenarios (Dolan and Mishra, 2013). Researchers are convinced that the age-dependent aerenchyma formation may not be the same for different root elongation rates, even though the patterns of distance-dependent aerenchyma formation are exactly the same; and the W₀-Gompertz model will be useful for describing longitudinal patterns (Yamauchi et al., 2021).

However, regardless of whether the elongation rates among different roots were the same, we could still not conclude that the aerenchyma formation had the same pattern at the same position from the root apex or at the same time points after root cell emergence within a root system. In this study, we observed that aerenchyma formation varied at the same distance from the root apex in rice roots with different lengths; and the longer the root was, the higher the aerenchyma percentage (Figure 1 and Supplementary Figure 1). The maximum aerenchyma percentage did not increase if the root length was >13 cm on the plateau (Supplementary Figure 1). Therefore, the longitudinal pattern of aerenchyma formation in rice adventitious roots must be closely related to the root length, and it is practicable to unify the convinced pattern by mathematically modeling in the development process of adventitious roots with different root lengths.

Mostly, the selection of non-linear models has been performed using the goodness-of-fit measures such as R², AIC, and BIC (Lucio et al., 2015; Bem et al., 2017), but these measures do not assess the biases of the parameters and may lead to selecting models that do not describe the biological phenomena correctly (Bates and Watts, 2007; Sari et al., 2019). Actually, when fitting the longitudinal pattern of aerenchyma formation with four types of non-linear models, all models had a good fit with an R² > 0.99 (Supplementary Table 1). However, the lower difference between the average mean of the bootstrap parameter estimates and the OLS estimates and the symmetrical distribution of these estimates were only observed in the Ti-Gompertz model rather than the W₀-Gompertz and other tested models (Figures 3, 4). Thus, the Ti-Gompertz model was the best for describing the longitudinal pattern of aerenchyma formation with different root lengths. Finally, according to the linear correlation of the parameters with the root length, we proposed the normalized Ti-Gompertz model as the following equation:

where AR is the percentage of the aerenchyma area to the cross-sectional area in the rice roots, D is the distance from the root apex, and L is the root length.

With the longitudinal patterns displayed by different root lengths, it was generally observed that longer roots displayed more developed aerenchyma than shorter roots in Eq. (11) and Supplementary Figure 2.

Parameter estimation can be defined as “a discipline that provides tools for the efficient use of the data in the estimation of constants appearing in mathematical models and for aiding in modeling of phenomena” (Dolan and Mishra, 2013). However, the two Gompertz forms were essentially the same model in mathematics (Tjørve and Tjørve, 2017). Therefore, we switched the equation form from the W₀-Gompertz model to the Ti-Gompertz model and reparameterized the Ti-Gompertz model. Apparently, compared with the W₀-Gompertz model, the Ti-Gompertz model provided parameter estimates close to being unbiased and good linear approximation; hence, this conversion made it possible to correlate root aerenchyma formation with any length of adventitious roots (Eq. (11); Figure 2).

Furthermore, determining the critical points and their biological meanings of models (e.g., Gompertz and logistic) is important in agricultural work (Ersoy et al., 2006; Prado et al., 2013; Pinho et al., 2014; Archontoulis and Miguez, 2015). The value at the IP of the sigmoid curve of the Gompertz model is locked at 36.8% of the upper asymptote (MVP) (Pinho et al., 2014; Tjørve and Tjørve, 2017). Around the IP, from the MAP to the MDP, the aerenchyma developed rapidly (Figures 2, 6). The asymptote increased and the IP decreased as the root length increased, also demonstrating that the development of aerenchyma in longer roots was stronger than that in shorter roots at the same position from the apex (Figure 2 and Supplementary Figure 1).

Several studies reported that the development of root aerenchyma could reduce root metabolic costs, and roots with large amounts of aerenchyma had lower respiration rates on a volume basis (Colmer, 2003; Zhu et al., 2010; Postma and Lynch, 2011; Saengwilai et al., 2014; Chimungu et al., 2015). Therefore, the aerenchyma formation improved plant performance under suboptimal nutrient and water availability and had a positive effect on overall root growth (Evans, 2004; Zhu et al., 2010; Jaramillo et al., 2013; Schneider et al., 2017). We described five stages of the Ti-Gompertz model based on the critical points and found that the speed of Stage III in the development of aerenchyma reached a maximum of 5.3–15.7% per cm and remained at 1.38–1.76 cm (Figure 6 and Table 1). It can be concluded that the cortical cells tended to die rapidly for the aerenchyma formation rather than for the respiration cost during this stage, and many anabolism processes were active. Thus, the cortical aerenchyma formation of rice adventitious roots at different longitudinal positions rather than that of the entire root might be a useful and precise trait for physiological functions such as gas transportation and water and nutrient acquisition. Eq. (11) of the Ti-Gompertz model, which we proposed in this study, helped to deeply understand the relationship between the trait and function of anatomy for rice adventitious roots.

To pinpoint the usefulness of the proposed model and the given results, we focused on some applications to the real data, and the volume of the entire aerenchyma could be calculated using the integral function of the Ti-Gompertz model (Figure 7). The method reported previously for the volume of the intercellular gas-filled spaces in roots measured the root porosity, which can be further enhanced by the aerenchyma formation (Jensen et al., 1969; Colmer, 2003). This study also reported alternative mathematical modeling for conveniently estimating the volume of the aerenchyma of the entire roots.