Root architecture characterization in relation to biomass allocation and biological nitrogen fixation in a collection of European soybean genotypes

This is anOpe Abstract – Soybean [Glycine max (L.) Merr] is the legume with the largest cultivated area worldwide and its yield depends largely on symbiotic nitrogen fixation and root architecture. This study aimed to explore the genetic variability of root architectural traits and di-nitrogen fixing activity in a small collection of nine European cultivars belonging to the same maturity group during their early stages. New image analysis approaches were implemented to characterise root architecture at high throughput. Significant genetic variability was identified for the width of the root system, root density, and for nitrogen fixation. This study allowed us to highlight trade-offs among root and nodule traits, and structural and functional traits. Finally, both the image analysis approach and the results could be used for breeding programs of soybean, that could take into account the root system architecture, when the plant interacts in symbiosis with N2-fixing bacteria.


Introduction
Soybean (Glycine max) is one of the most cultivated oilseed crop worldwide. Despite the European and French policy to achieve protein autonomy, the production of soybean (2 813 260 tonnes between 2019 and 2020), is still far away from the level of the importation (33 264 815 tonnes during the same period, European Commission's agricultural and rural development department, 2021). These imports have accelerated with the intensification of farming activities (de Visser et al., 2014) and, to a lesser extent, with the transformation of oils into biogas (Westhoff, 2009). Europe is thus striving to intensify its soybean production, in particular because of agroecological advantages related to this crop. First, this leguminous crop does not require nitrogen fertilizer thanks to biological nitrogen fixation, and on the contrary enriches the soil in nitrogen with its nitrogen-enriched residues (Peoples et al., 2009;Hungria and Mendes, 2015). This culture requires little or no phytosanitary treatment (Reijnders and Soret, 2003) and produces seeds which are rich in oil, protein, fibre, starch and other essential nutrients making it valuable for the production of food and feed (Preissel et al., 2015;FoodData Central). However, soybean yields per unit area have very little increased during the last 10 years, suggesting that more research efforts for soybean breeding should be engaged , especially in the current context of climate change.
Over the last decades, many studies focused on the mechanisms regulating the mutually beneficial association between soybean and N 2 -fixing rhizobacteria (Bradyrhizobium japonicum), that is involved in biological nitrogen fixation (BNF). While bacteria provide a source of nitrogen to the plant, plants provide in return carbon sources for nodule growth and for rhizobial metabolism. In contrary to temperate legumes for which the nodule development is indeterminate, in soybean, BNF occurs within nodules of determinate growth whose shape is spherical (Popp and Ott, 2011;Udvardi and Poole, 2013).
More recently, in order to identify new target traits for breeding, especially for drought tolerant purposes, several studies tackled the genetic variability of the root system architecture in soybean. Compared to other agronomically traits of interest, the selection of new soybean varieties based on root architecture criteria has lagged (York et al., 2015). This absence of selection towards an optimal root system, for example, has led to an impoverishment of the root diversity of American varieties (Falk et al., 2020). However, it is now known that a deep root system is beneficial for the uptake of nitrogen and water in deep layers and that, in contrast, a shallow system with many adventitious root is beneficial for the uptake of relatively immobile nutrients such as phosphorus (Ho et al., 2005;Lynch, 2011;He et al., 2017).
The genetic variability of root architecture in soybeans has been evaluated in different maturity groups, and has highlighted a significant variability of root traits during their early growth (Falk et al., 2020;Dayoub et al., 2021). Some of these variables have also been correlated with aerial traits, such as the root length which was positively correlated with shoot dry matter (Dayoub et al., 2021). These studies provided very interesting results by revealing the root growth potential of each of the genotypes. However, the influence of symbiosis related to BNF on root system architecture in legumes is still poorly understood (Concha and Doerner, 2020), while it seems necessary to be able to quantify the extent of the trade-off between root growth and nodular growth. Indeed, the competition for carbon between roots and nodules, can limit the growth and development of any of the organs, and even their functioning (Voisin et al., 2003). This was sustained by the work from Yang et al. (2017) that highlighted several root traits were significantly affected by the presence or absence of rhizobia in the soil.
In this study, we therefore wanted (i) to assess the nodulated root architecture of several French soybean varieties, when the acquisition of N was solely based on BNF, and (ii) to characterize the efficiency of BNF during their early growth.
Seeds were calibrated and pre-germinated for 4 days at 20°C in a Fitoclima S600 germinator (Aralab, Rio de Mouro, Portugal) before transplantation. Seedlings were grown in RhizoTubes ® filled with a mixture of perlite and sand (3:1, v:v), allowing the visualization of the root system (Jeudy et al., 2016). Two seedlings were then transplanted in each of the 38 RhizoTubes ® , consisting in 8 to 12 RhizoTubes ® per genotype. At the same time, each seedling was inoculated with 1 mL of Bradyrhizobium diazoefficiens MIAE00426 (MSDG1996, G49; (Lagacherie et al., 1977)) obtained from the collection "Microorganismes d'intérêt Agro-Environnemental" located at UMR Agroécologie, INRAE, Dijon, France and corresponding to 10 8 rhizobia. A second inoculation following the same characteristics of the first one was performed two days after the transplantation. Plants were grown in a greenhouse of the Plant Phenotyping Platform for Plant and Microorganism Interaction (4PMI) at INRAE in Dijon (France) (47°32 0 N, 5°02 0 E) during 23 days. Mean day/night temperatures were 18.5/15°C, the relative humidity was controlled to 45% and the photoperiod was set to 16 h thanks to an artificial lighting (PAR of 272 mmol.m À2 .s À1 ) supplied with sodium lamps (400 W lamp, HPS Plantastar, OSRAM, Munich, Germany). Plants were watered 3 times a day with 250 ml of N-free nutritive solution as described in Voisin et al. (2003).
Each RhizoTube ® was imaged three times a week using the high throughput aerial and root phenotyping booths in the RhizoCab ® HD (Jeudy et al., 2016).

Image analysis
The 342 high-resolution images of the growth kinetics of soybean root development were analysed by combining two methods. First, high-resolution images (12 000 * 12 000 pixels) ( Fig. 1) were divided into 16 smaller images (3000 * 3000 pixels). Then, on each small image, a segmentation was performed with the machine learning software Ilastik (Berg et al., 2019). The small-segmented images were then reassembled and analysed using a program written under python programming language. This program selected the plant to be analysed (2 plants per image), transformed the image into a matrix and counted the number of black and white pixels to measure the root-projected area (Fig. 1). The submodule dedicated to morphology measurements, from "skimage" package, as well as the sub-module patches of the "matplotlib" package allowed to measure the area and the perimeter of the convex shape around the root system as well as the width and the height of this convex shape corresponding to the width and the depth of the root system (Fig. 1). The root density was calculated as the ratio of the root projected area to the convex shape area. Python scripts used to perform these analyses can be found in the Supplementary Material Appendix A1.
Using the Fiji software's multi-point function (Schindelin et al., 2012), the coordinates of the orthonormal frame (x, y) of the image of the rightmost point, the leftmost point of the root system as well as the branching on the taproot and the tip of the taproot were measured during the growth of the plant for each plant. Then, the distance y between the two outermost points was calculated and corresponded to the width of the root system, and the distance x between the first branch and the tip of the taproot corresponded to the depth.
The characterization of nodule traits was performed directly on the high-resolution images, using the Cell Counter plugin of the Fiji software. This plugin allowed to categorize After segmentation, each pixel of the image was worth the value 0, 1 or 2 and corresponded respectively to the background, the roots or the nodules (black image). (D) Image conversion into white, black and red images (respectively the background, the roots and the nodules). These segmented images were then reassembled. (E) The two plants were extracted from the image to be analyzed individually. (F) The root projected area, nodule projected area were then measured as well as the width, depth, area and perimeter of the convex hull. Once the template is created it takes approximately 14 minutes from C2 to F. and count the nodules. Each nodule was marked with a particular-coloured point, in order to discriminate pink efficient nodules from the white inefficient ones, and the position of each nodule on the root system (insertion on the taproot or on lateral roots).

Plant measurements at harvest
After 23 days of growth in the RhizoTubes ® , all soybean plants were at the V1 stage and plants were harvested individually. The nodulated roots were separated from the aerial parts. Nodules were then manually removed from the root system. The number of lateral roots inserted on the taproot was counted. Each plant compartment (roots, nodules and shoots) was dried for 48 h at 80°C before dry biomass (BM) measurement. Then, plant tissues were individually ground into a fine powder (ball mill Retsch MM 400) for carbon and nitrogen concentration analyses, using an elemental analyser (Thermo electron NC2500, Courtaboeuf, France).
The total amount of dry biomass (BMtotal) was calculated as the sum of the dry biomasses of all plant compartments (roots, nodules, shoots). Then, in order to observe the allocation of dry biomass within the plant, shoot biomass to below-ground biomass ratio and nodule biomass to nodulated root biomass ratio were calculated.
The amounts of C and N contained in the whole plant (QC and QN respectively) were calculated following: [N] i and BM i are C and N concentrations and dry biomass, respectively, in each plant compartment i.
The plant C:N ratio was calculated as: Nodule N 2 fixation specific activity, expressed in gN. gBM nodule À1 .d À1 was calculated following the method described in Prudent et al. (2016).

Statistical analyses
All statistical analyses were performed with R software (R Core Team, 2019). The evolution of the number of nodules over plant growth was fitted to a three parameters Gompertz curve (Pegelow et al., 1977;Tjørve and Tjørve, 2017) following this equation: where A is the asymptote corresponding to the maximum number of nodules, Kip is the nodulation rate at the inflexion point and Tip is the time at the inflection point of the curve.
The Gompertz parameters A, Kip and Tip were estimated using the nls function in R software.
For each model of each genotype, the goodness-of-fit was evaluated with the relative root mean squared error (RRMSE) criterion (Kobayashi and Salam, 2000): where O i is the observed value, S i is the corresponding predicted value, N is the number of observed data, and O mean is the mean of all measured values. The smaller the value of RRMSE, the better the goodness of fit. For each measured or calculated trait, a global comparison test was conducted using the non-parametric Kruskal's test followed by Tukey's HSD (Honest Significant Difference) test. For the traits corresponding to ratio, they were first log10 transformed before statistical analysis. For the three parameters of the Gompertz curve, a bootstrap sampling method was used before statistical analysis.
Principal component analysis (PCA) was performed (16 traits) using the "FactoMineR" package and graphed using the "factoextra" package in R software. Intergenotypic Pearson's correlations were calculated and considered as significant at the 0.05 probability level.

Results and discussion
In order to quantify the genetic variability of root and nodule traits, we decided (i) to analyse structural (including biomass allocation, morphological characterization) and functional (N 2 fixation) traits and (ii) to focus on a small collection of nine varieties belonging to the same maturity group. This choice was motivated by the need to compare root traits from plants at the same developmental stage, and was supported by a first study conducted by Dayoub et al. (2021). The full dataset is provided in Supplemental Material Table S1.
3.1 Genetic variability of biomass allocation was observed among varieties belonging to the same maturity group First, the whole plant dry biomasses were compared among genotypes when they reached the V1 stage (Fig. 2). They ranged from 362 mg in Straviata to 505 mg and 516 mg in Wendy and Sinema, respectively. Interestingly, for Sinema, its higher biomass was partly due to a higher nodule biomass than that of the other genotypes (65% higher than the average of the other genotypes). In order to get a better insight of the allocation of biomass between the three main organs (shoots, roots and nodules), ratio of biomasses were calculated, and revealed that the shoot to nodulated root ratio (Fig. 3A) was similar among genotypes, except for Isidor. Isidor allocated 1.5 more biomass to shoots than the eight other genotypes. This particular feature of Isidor was related to a higher shoot length (data not shown), that could be an advantage for light competition with weeds (Colbach et al., 2020), but the low amount of biomass dedicated to nodulated roots could hamper soil mineral resource uptake under limiting conditions. When focusing on the nodulated root system, we highlighted that the relative allocation of biomass to nodules within the root system significantly differed among genotypes (Fig. 3B). The two extreme genotypes were Sinema and Wendy: for a same quantity of below-ground biomass, Sinema allocated 1.75 times more biomass to nodules than Wendy. Under limiting soil resource conditions (drought, nutrient deficiency...), it is tempting to speculate that Sinema could favor nitrogen acquisition via BNF at the expense of water or soil nutrient uptake.

Root architectural traits and N 2 -fixation activity
The use of the innovative device RhizoTubes ® (Jeudy et al., 2016) associated with specific image analysis (Fig. 1) allowed us to easily characterize root system architecture of the nine varieties over time. These root architectures were highly contrasted, depending on the genotype as shown in Figure 4. Although the length of the taproot was not a discriminant trait, the width of the root system together with the root density remarkably discriminated genotypes (Fig. 5). In contrast to our study, some authors highlighted a genetic variability for root depth at early stages (Manavalan et al., 2010;Matsuo et al., 2013), or during flowering (Fried et al., 2018), but considered genotypes from contrasted maturity groups. In our study, Sumatra genotype displayed a wider root system compared to Wendy (24.6 cm and 15.3 cm, respectively). Furthermore, for this trait, differences among genotypes were observable as soon as 16 days after transplantation (Fig. 5A). The more the plant developed, the more marked these differences were.
Regarding root density, as soon as 6 days after transplantation some genotypes displayed contrasted behaviours (Fig. 5C): Wendy (41%), Pallador (42%), Stocata (41%) and Steara (41%) had a higher root density than Sumatra (33%), Isidor (32%) and Sinema (30%), but this trend was not maintained until the end of the experiment. For all genotypes except Wendy and Sinema, a decrease in root density between 6 and 10 days was observed. These results underline that some traits (such as root density) cannot be measured very early during root development as they follow genotype-dependent kinetics, contrarily to other traits that seem to be established at a very early stage (width of the root system). Finally, the most contrasted genotypes in terms of root density were Wendy (52%) and Straviata (32%). Wendy could be a good candidate for soils containing low phosphorus content. Indeed, dense root systems are more adapted for this type of deficiency (Al-Ghazi et al., 2003;Lynch, 2011;Lynch, 2019).
The nodulation process was also quantified by following the number of nodules over time for each genotype (Fig. 6A). Although the dry biomass allocated to nodules differed according to the genotype (Fig. 2), the number of nodules was similar whatever the genotype at the end of the experiment (Fig. 6B). By fitting the number of nodules over time to a threeparameters Gompertz model, we were able to statistically distinguish the genotypes according to the different parameters of the model (Tab. 1). The Kip parameter, which reflects the maximal rate of nodule development, was lower for Steara, than for Sinema for which the value of the parameter doubled. Concerning the parameter Tip, which represents the time when Kip was the highest, it highlighted that some genotypes initiated earlier their nodulation than others (12 days for Sinema and Isidor, 16 days for Sinema), and could thus be Fig. 2. Dry biomass of the nine soybean genotypes, measured at harvest (23 days after transplantation). Letters indicate statistically significant differences among genotypes for plant total dry biomass (blue letter) or for each organ (within the column) (Tukey'ś HSD test. P < 0.05). more responsive to plant N needs. In pea, this behaviour has been shown to be an advantage to allow a better plant recovery following a drought period (Couchoud et al., 2020).
Genotypes also differed regarding the location of their nodules: while more than 80% of the nodules were located on the taproot (primary root) in Sinfonia, they were only 40% in Steara. To our knowledge this feature has never been attributed to an ability to tolerate specific edaphic conditions, but we can speculate that a distribution of nodules throughout the root system could help to face heterogenous soil constraints (e.g. anoxic conditions...), unlike a limited nodulation zone.
The specific activity of nodules was also quantified and revealed strong differences among genotypes (Fig. 6E). This N 2 -fixing activity can be related on one hand to the efficiency of the rhizobial strain (Senaratne et al., 1987) and on the other hand to the interaction between the plant and the rhizobia, leading to more or less substrate for nitrogenase activity. In our study, the same rhizobial strain was used for all the genotypes, meaning that the observed differences only originated from the soybean genotype. Sumatra was the less efficient for BNF, while Isidor has a three times higher N 2 -fixation activity, followed by Pallador and Wendy. Knowing BNF is highly dependent from abiotic stresses such as drought (Streeter, 2003;Fenta et al., 2014;Sharaf et al., 2019), temperature (Munevar and Wollum, 1982;Keerio et al., 2001), salt (Singleton and Bohlool, 1984), soil acidity (Hungria and Vargas, 2000), or nutrient deficiency like iron (Brear et al., 2013), the plasticity of BNF to some of these conditions could be evaluated for the genotypes showing the most contrasting specific activities (i.e. Isidor, Pallador and Wendy).

Correlations among root and nodule traits
A principal component analysis was performed in order to summarize differences observed among genotypes (Fig. 7). Traits that accounted the most to axis 1 were the coefficients of the Gompertz curve A and Tip corresponding to the maximum number of nodules predicted by the model, and the time when the increase in nodules is the highest, respectively, and the width of the root system. For axis 2, traits were root density, percentage of nodules located on the primary root and Kip, the nodulation rate of the number of nodules. Based on our nodulated root characterization, three main groups of genotypes could be identified: a first one composed by Wendy, Sinema, Pallador, Isidor and Sinfonia was characterized by a high root density, and a high rate of nodulation; a second group composed by genotypes Straviata, Steara and Sumatra was mainly characterized by a growth in width of the root system and a third group composed by genotype Stocata was characterized by an average behaviour between the two first groups.
Intergenotypic correlation analysis among traits was performed in order to better understand the trade-off between roots and nodules (Fig. 8). We observed that total dry biomass of the plant correlated with several traits related to root and nodules, as observed by Yang et al. (2017) but the determinism of the nodulation process still needs to be deciphered. In pea, the nodulation process was shown to be tightly adjusted to plant growth rate, in order to fulfil nitrogen needs of the plant (Voisin et al., 2010). Then, nodule specific activity was negatively correlated to the relative allocation of dry biomass to nodules within the root system, suggesting that genotypes that invested in structural component of N 2fixation (nodule dry biomass) allocated less energy for nitrogenase activity. From an architectural point of view, the number of lateral roots was negatively correlated to the rate of nodulation illustrating the trade-off between root branching and nodule initiation, as already suggested by pioneer work by Nutman (1948).  Root system architecture over time in the nine soybean genotypes: Root system width between 6 and 23 days after transplantation (A) and at harvest (B). Root density from 6 to 23 days after transplantation (C) and at harvest (D). The 95% confidence interval with local regression adjustment was drawn in transparency. In the boxplots, letters indicate significant differences among genotypes (Tukey'ś HSD test. P < 0.05).  Number of nodules from 6 to 23 days after transplantation (A) and at harvest (B). Percentage of nodules located on the primary root from 6 to 23 days after transplantation (C) and at harvest (D). Mean specific activity of nodules during the 23 first days of growth (E). The 95% confidence interval with local regression adjustment was drawn in transparency. In the boxplots, letters indicate significant differences among genotypes (Tukey'ś HSD test. P < 0.05).

Conclusion
Until now, soybean varieties have been classified into different maturity groups depending on their needs in terms of heat and photoperiod: thirteen maturity groups exist throughout the world (Poehlman, 1987). In the present study, the analysis of root and nodule traits allowed us to identify a genetic variability associated on one hand to the acquisition of nitrogen from BNF, and on the other hand to the potential water and nutrient uptake, even within a same soybean maturity group. Three groups of genotypes were identified according to their contrasted morphology and nitrogen fixing activity. We were also able to highlight trade-offs among root and nodule traits, and structural and functional traits. But nodulated root architecture of soybeans depends also on the environment. Indeed, the architecture can be modified by nutrient concentration (Linkohr et al., 2002;Li et al., 2016), soil structure, water availability (Hoogenboom et al., 1987;Fenta et al., 2014;Xiong et al., 2020), pollutants (Bashir et al., 2019)... defining their plasticity. Now, further studies should consider the relationship between root architectural traits, nodulation and water and nutrient uptake under various soil conditions. This should include the influence of nodulated root traits on soil microbial communities (Hetrick, 1991;Saleem et al., 2018) that play a critical role in bio-geochemical cycles, and thus on nutrient availability for plant.

Supplementary Material
Appendix 1. Python scripts used for image analysis. Table S1. Dataset summarizing measured and calculated traits on the nine soybean genotypes.
Acknowledgements. This work was supported by the ECODIV project of the Plant2Pro Carnot Institute. The authors are grateful to Cécile Révellin and Baptiste Serbourse for rhizobial preparation, to Vincent Durey, Sylvie Girodet and Céline Latour for their technical support during harvests, and to the members of the High Throughput Phenotyping Platform the 4PMI regarding their support during plant growth. The authors also thank the breeding companies that supplied the soybean seeds. Fig. 7. Principal component analysis of the nine soybean genotypes. The two principal components accounted for 27.5% (PCA1) and 16.6% (PCA2) of the total variation. Ellipses in color represent the genotypes and the color of the arrows indicates the contribution of each trait: a color gradient was used from red (high contribution) to blue (low contribution). Additional variables that were not considered to create the PCA were colored in blue.