<figure><img src="image_1.png" /></figure>

HAL
open science

# Genetic control of rhizosheath formation in pearl millet

Carla de la Fuente Cantó, Mbourou Nelly Diouf, Papa Mamadou Sitor Ndour, Marilyne Debieu, Alexandre Grondin, Stephanie Passot, Antony Champion, Carlos Barrachina, Marine Pratlong, Pierre Gantet, et al.

<figure><img src="image_2.png" /></figure>

## To cite this version:

Carla de la Fuente Cantó, Mbourou Nelly Diouf, Papa Mamadou Sitor Ndour, Marilyne Debieu, Alexandre Grondin, et al. Genetic control of rhizosheath formation in pearl millet. Scientific Reports, 2022, 12 (1), pp.9205. 10.1038/s41598-022-13234-w. hal-03714569

HAL Id: hal-03714569

https://hal.inrae.fr/hal-03714569

Submitted on 5 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

Check for updates

OPEN

# Genetic control of rhizosheath formation in pearl millet

C. de la Fuente Cantó<fnref n="1" />, M. N. Diouf²,³,⁷, P. M. S. Ndour²,³, M. Debieu<fnref n="1" />, A. Grondin<fnref n="1" />,⁴,⁵, S. Passot<fnref n="1" />, A. Champion<fnref n="1" />, C. Barrachina⁶, M. Pratlong⁶, P. Gantet<fnref n="1" />, K. Assigbetsé²,³, N. Kane⁴, P. Cubry<fnref n="1" />, A. G. Diedhiou⁴,⁷, T. Heulin⁸, W. Achouak⁸, Y. Vigouroux<fnref n="1" />, L. Cournac² & L. Laplaze<fnref n="1" />,⁴✉

The rhizosheath, the layer of soil that adheres strongly to roots, influences water and nutrients acquisition. Pearl millet is a cereal crop that plays a major role for food security in arid regions of sub-Saharan Africa and India. We previously showed that root-adhering soil mass is a heritable trait in pearl millet and that it correlates with changes in rhizosphere microbiota structure and functions. Here, we studied the correlation between root-adhering soil mass and root hair development, root architecture, and symbiosis with arbuscular mycorrhizal fungi and we analysed the genetic control of this trait using genome wide association (GWAS) combined with bulk segregant analysis and gene expression studies. Root-adhering soil mass was weakly correlated only to root hairs traits in pearl millet. Twelve QTLs for rhizosheath formation were identified by GWAS. Bulk segregant analysis on a biparental population validated five of these QTLs. Combining genetics with a comparison of global gene expression in the root tip of contrasted inbred lines revealed candidate genes that might control rhizosheath formation in pearl millet. Our study indicates that rhizosheath formation is under complex genetic control in pearl millet and suggests that it is mainly regulated by root exudation.

**Abbreviation**

RAS Root-adhering soil

Pearl millet is a small-seeded tropical cereal that was domesticated about 4500 years ago in the Sahelian part of West Africa<fnref n="1" />. It is mostly grown in dry and poor soils as a rainfed crop and is therefore well adapted to environments prone to drought and heat stress for which it harbours largely untapped genetic diversity in the locally adapted cultivated and wild pearl millets<fnref n="1" />⁻³. The outstanding capacity for growing in harsh environments highlights the great potential of pearl millet as a biological model to investigate crop adaptation and resilience to abiotic constraints, as well as its key role for food security in some semi-arid tropical regions in Africa and Asia. Still, pearl millet yield remains low for two main reasons: the difficulty to reach its potential yield in constrained environments and the little attention that the crop has received from breeding programmes³.

Root traits are emerging as new targets for breeding more sustainable and resilient crop varieties in global climate change scenarios⁴. The root system is responsible for plant water and nutrient acquisition. Phenotypic selection of root ideotypes combining architectural, anatomical and physiological traits has been proposed as a way to optimise access to soil resources in specific agroecosystems and crop management practices⁴. Besides root architecture, anatomy and physiology, the rhizosphere, the volume of soil around the root influenced by the root⁵, can be regarded as a plant extended phenotype and therefore a target for breeding more sustainable crops⁶,⁷. Indeed, the dynamic interplay between root, soil and microbiota in the rhizosphere eases adaptation to changing environments and can have a remarkable impact on plant fitness⁷⁻⁹. The intricate relationships in the rhizosphere define a belowground niche where soil moisture, organic matter content, the composition of the microbial community and its activity are different from the bulk soil¹⁰,¹¹. Plants benefit from this interaction especially in constrained environments where access to nutrients and water is restricted¹².

<footnote n="1">¹DIADE, Université de Montpellier, IRD, CIRAD, Montpellier, France. ²Eco&Sols, Université de Montpellier, IRD, CIRAD, INRAE, Institut Agro, Montpellier, France. ³Laboratoire Mixte International Intensification Écologique Des Sols Cultivés en Afrique de L'Ouest (IESOL), Dakar, Senegal. ⁴Laboratoire Mixte International Adaptation Des Plantes Et Microorganismes Associés Aux Stress Environnementaux (LAPSE), Centre de recherche de Bel-Air, Dakar, Senegal. ⁵CERAAS, Thiès, Senegal. ⁶Montpellier GenomiX, Montpellier, France. ⁷Département de Biologie Végétale, Faculté Des Sciences Et Techniques, Université Cheikh Anta Diop, Dakar, Sénégal. ⁸Aix Marseille Univ, CEA, CNRS, BIAM, LEMIRE, Laboratory of Microbial Ecology of the Rhizosphere, ECCOREV FR 3098, 13108 Saint Paul-Lez-Durance, France. ✉email: laurent.laplaze@ird.fr</footnote>
The rhizosheath size, or root-adhering soil mass, is a proxy in the study of this complex extended phenotype and an interesting potential target for breeding programmes¹³. Rhizosheath formation was first noticed as the sandy sheath surrounding the roots of desert plants¹⁴ and later reported to occur across many angiosperm orders¹⁵. Increased rhizosheath size has been correlated with enhanced wheat and foxtail millet performance in drying soils¹⁶,¹⁷. In barley and oat, rhizosheath formation has been related with improved acquisition of major and essential trace elements in limiting water conditions¹⁸,¹⁹. A combination of root architectural and anatomical traits and the secretion of root exudates and mucilage have been connected to soil aggregation to the root¹³,²⁰. For instance, root branching, root hair formation or symbiosis with arbuscular mycorrhizal fungi (AMF) have been associated to some extent with rhizosheath establishment¹⁵,¹⁷,²¹. Root architectural traits have been found crucial for rhizosheath formation in wheat and foxtail millet¹⁷,²². On the other hand, root exudates composition and mucilaginous polymers released by root-associated microorganisms impact the stability of soil aggregates that bind around the root²³. Root growth and exudates exert dynamic changes in the rhizosphere physical properties and hydraulic processes that affect soil nutrient dynamics and the composition of the rhizosphere associated microbiota⁹,²⁴⁻²⁶. Despite the inherent complexity linked to the effect of exudates in the rhizosphere, some studies showed their direct relationship with rhizosheath formation. For example, greater mass of mucilage exuded by chickpea roots were linked with the formation of larger and more porous rhizosheaths capable of storing more soil moisture in drought tolerant cultivars²⁷. In annual crops such as wheat, barley and maize, there is evidence of the remarkable plant genetic influence in the formation of rhizosheath and the processes of rhizodeposition influencing rhizosphere microbial activities¹⁸,²⁸⁻³⁰, however few studies have aimed to dissect the genetics underlying the conformation of this extended root phenotype¹⁸,²⁸,³¹,³².

In previous studies, we reported a remarkable genotypic variability for root-adhering soil aggregation in pearl millet³³. Moreover, this variability was associated with changes in rhizosphere microbiota structure and function³³,³⁴. Here, we analysed the relative contribution of root architectural characteristics and root colonization by AMF on root-adhering soil aggregation in pearl millet. We then combined a genome wide association analysis (GWAS), with bulk segregant analysis (BSA) and transcriptomic data to dissect the genetic bases of this complex trait.

## Results

**Root-adhering soil aggregation is weakly correlated to root hair traits in pearl millet.** Several root traits have been proposed to contribute to root-adhering soil aggregation (as an integrative phenotype) including root hair development, root architecture, and arbuscular mycorrhizal symbiosis. We therefore analysed the contribution of these different traits to root-adhering soil aggregation in pearl millet. For this, we analysed correlation between root-adhering soil aggregation, root architecture, root hair length and density and frequency and intensity of root colonization by AMF in eight inbred lines with contrasted rhizosphere aggregation phenotype³³ after four weeks of growth (Fig. 1A). Among root architecture traits, only the average root diameter (AvgDiam) was weakly and negatively correlated ($p=0.012, r^2=0.057$) with root-adhering soil aggregation (Fig. 1B, Table 1). For root hairs, only average length of root hairs (AvgLRH) was weakly and positively correlated to root-adhering soil aggregation ($p=0.005, r^2=0.077$; Fig. 1C, Table 1). No significant correlation was observed for all other traits including frequency and intensity of root colonization by AMF (Fig. 1A, Table 1). Similar results were found in two independent experiments (2018 and 2020; Supporting Information Table S1).

Altogether, our results suggest that root hairs development could play a weak role in root-adhering soil aggregation in pearl millet and that root architectural traits and AMF colonization rate have no significant impact.

**Genetic bases of rhizosheath formation in pearl millet.** We previously reported the phenotyping of a panel of pearl millet inbred lines for root-adhering soil aggregation³³. Briefly, a total of 1408 plants corresponding to 181 inbred lines were phenotyped and we recorded an almost four-fold variation in rhizosheath size (RAS/RT ratio), ranging from 7.4 (ICML-IS 11139) to 26.3 (ICML-IS 11084)³³. Here, we used these data to evaluate the heritability of root-adhering soil aggregation. A broad sense heritability of 0.72 was computed, suggesting that root-adhering soil aggregation is largely under genetic control. Altogether, these data indicate that root-adhering soil formation has high heritability and that a large genetic diversity exists in pearl millet.

We therefore analysed the genetic bases of root-adhering soil formation using association genetics. Out of the 181 inbred lines, 139 lines with good quality data for phenotype and genotype were retained to perform the association study. As a first step, we conducted a population structure analysis of the 139 lines (Supplementary Fig. S1) that confirmed the negligible genetic structure previously reported for this panel³⁵. A total of 381,899 SNPs was used for the association analysis. We first calculated the least square means of the trait root-adhering soil aggregation (RAS/RT ratio) across the different experiments. The ratio ranged from 12.4 to 26.3 with an average of 18.0. The LFMM model for GWAS identified 53 significant SNPs ($p$ value < 0.0001) across the genome (Fig. 2A), defining 34 significant regions or QTLs considering windows of 50 kb up and downstream of significant positions to define significant regions. The proportion of phenotypic variance accounted for the most significant SNPs defining QTLs ranged from 9.2 to 15.6% indicating that the corresponding QTLs had small phenotypic effect.

We compared these results with two other GWAS methods (Fig. 2B,C). Thirty-nine of these SNPs included in 25 QTLs defined through LFMM were also found significant with EMMA, and 19 significant SNPs assigned to 14 QTLs using MLM model in GAPIT (Supplementary Table S2). Fifteen SNPs in 12 QTL regions were found significant across the three GWAS methods. Altogether, our GWAS analysis revealed 12 potential QTLs controlling root-adhering soil aggregation in pearl millet.

To back up our GWAS analysis, we performed bulk segregant analysis in a F2 population derived from a cross between two lines with contrasted RAS/RT phenotypes, ICML-IS 11139 (low RAS/RT) and ICML-IS 11084 (high
<figure><img src="image_5.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_7.png" /></figure>

**Figure 1.** Relation between root soil aggregation, root architecture, root hair development and arbuscular mycorrhizal symbioses. (A) Pearson correlation between traits using adjusted lsmeans across two experiments conducted in different years. (B) Linear regression between root diameter and root soil aggregation. Points represent the mean value of the traits for inbred lines across the two experiments. (C) Linear regression between root hair length and root soil aggregation. Points represent the mean value of the traits for inbred lines across the two experiments.

RAS/RT; Fig. 3A). F2 plants were phenotyped in five consecutive blocks together with the parental inbred lines. We confirmed the contrasted RAS/RT ratio of the parental lines with average values of 15.0 (ICML-IS 11139) and 32.7 (ICML-IS 11084, Fig. 3B). Ten individual F2 lines were dropped from the analysis leaving 547 F2 with RAS/RT ratio ranging from 1.6 to 54.8 and with an average value of 22.3. The phenotype distribution of the F2 was slightly skewed towards high values of RAS/RT ratio and showed a significant block effect (Supplementary Fig. S2). We used log transformation of RAS/RT ratio in our analysis of variance and selected lines with extreme values of the residual term for the bulks. The bulks consisted in two groups of 55 F2 lines each, with RAS/RT ratio average values of 11.0 for the low RAS/RT bulk and 38.2 for the high RAS/RT bulk (Supporting Information Table S3). A total of 223.6 Mbp reads were mapped to the target enriched regions and used for SNP calling.
<table><thead><tr><th rowspan="2">Trait</th><th colspan="8">Root architecture</th><th colspan="2">Root hairs</th><th colspan="2">AM symbiosis</th></tr><tr><th>RAS/RT</th><th>L</th><th>RSA</th><th>AvgDiam</th><th>LFR</th><th>LTR</th><th>SAFR</th><th>SATR</th><th>AvgLRH</th><th>AgDRH</th><th>F%</th><th>1%</th></tr></thead><tbody><tr><td>RAS/RT</td><td>1</td><td>0.017</td><td>-0.081</td><td><strong>-0.239</strong></td><td>0.053</td><td>-0.112</td><td>0.037</td><td>-0.110</td><td><strong>0.278</strong></td><td>0.069</td><td>0.283</td><td>-0.133</td></tr><tr><td>p valuc</td><td>0</td><td>0.863</td><td>0.397</td><td><strong>0.012</strong></td><td>0.579</td><td>0.242</td><td>0.697</td><td>0.250</td><td><strong>0.005</strong></td><td>0.488</td><td>0.463</td><td>0.744</td></tr></tbody></table>

**Table 1.** Correlation matrix of root soil aggregation and root parameters. Value from 2 independent experiments on contrasted pearl millet lines using Spearman's correlation test. Ratio (RAS/RT) between the mass of root-adhering soil (RAS) and root tissue biomass (RT), Total root length (L), Root Surface Area (RSA), Average Root Diameter (AvgDiam), Total Length of Fine Roots (LFR), Total Length of Thick Roots (LTR), Surface Area of Fine Roots (SAFR), Surface Area of Thick Roots (SATR), Average Length of Root Hairs (AvgLRH), Average Density of Root Hairs (AgDRH), Frequency of mycorrhization (F%), Intensity of mycorrhization (I%). Significant values are in bold.

<figure><img src="image_8.png" /></figure>

**Figure 2.** Genome-wide association studies (GWAS) for rhizosheath size in pearl millet. Manhattan plots and QQ plots obtained with three GWAS methods. (A) Efficient Mixed Model Association or EMMA, (B) Mixed linear model or MLM and (C) Latent Factor Mixed Model or LFMM. Each Manhattan plot shows the -log10 *p* value of the statistic (y axes) for each SNP position (x axes). The dashed line delimits the threshold for highly significant SNPs (*p* value < 10⁻⁴).

We identified a group of 23,160 SNP variants (1.5 SNPs per 100 kb in average) between the bulks. The average sequencing depth was high with 887X and 863X in the small and high RAS/RT bulk respectively. The NGS-based BSA analysis revealed significant differences in the allele frequency of 380 SNPs at the 95% confidence interval (Table 2, Fig. 3C). These SNPs defined five significant chromosome regions linked to the segregation of the RAS/RT ratio phenotype: three on chromosome 5 (RAS5.1, RAS5.2 and RAS5.3) and two on chromosome 6 (RAS6.1,
<figure><img src="image_9.png" /><figcaption>**Figure 3.** Genetic dissection of root soil aggregation in pearl millet by Bulk Segregant Analysis (BSA). (A) Cross established for Bulk Segregant Analysis (BSA) between two pearl millet inbred lines with contrasted rhizosheath phenotype. (B) Boxplot showing the distribution of RAS/RT ratio in line ICML-IS 11139 (N = 29), ICML-IS 11084 (N = 27) and F2 population (N = 547). (C) Comparison between GWAS and BSA results. Top figure represents the Manhattan plot of the GWAS by LFMM ridge method³⁶. The x-axis corresponds to the position of the 381,899 SNPs identified by GBS in a group of 139 inbred lines. The vertical axes correspond to the -log10 *p* value of the statistic. The dashed line delimits the threshold for highly significant SNPs (*p* value < 10⁻⁴). Bottom figure shows the significant regions associated with root soil aggregation identified by BSA using bulks of contrasted F2 lines from a bi-parental cross. The plot shows the Euclidean Distance statistic profile (y axis) across the seven pearl millet chromosomes (x axis). The dashed line indicates the 95% confidence interval threshold for the localisation of significant regions. In both plots, the shaded area delimits the extent of the five significant regions identified by BSA and the overlap with significant SNPs identified by GWAS and the correspondence with the BSA peaks found.</figcaption></figure>

RAS6.2; Table 2). The smallest genomic region defined corresponded to RAS5.3 with 10.6 Mbp and 41 significant SNPs. In contrast, the largest significant region corresponded to RAS5.2, with 45 Mbp and 307 significant SNPs.

With this information we looked for chromosome regions coming up significant across analysis. We compared the mapping position of QTL regions defined by the significant SNPs identified in the GWAS analysis and the significant regions identified using BSA in the F2 population. Interestingly, the range of four out of five BSA significant regions was found to overlay with the position of significant SNPs defined by GWAS (Fig. 3C,
<table><caption>**Table 2.** Significant genomic regions identified by Bulk Segregant Analysis (BSA) for root-adhering soil aggregation (i.e. RAS/RT) at the 95% confidence interval. ¹Position of the most significant SNP in the region range. ²Limits of the significant region considering the overlapping confidence interval of significant markers in the region.</caption><thead><tr><th></th><th>Chr</th><th>Peak position (Mbp)<sup>1</sup></th><th>Region range (Mbp)<sup>2</sup></th><th>Region length (Mbp)</th><th>Number sig SNPs</th></tr></thead><tbody><tr><td>RAS5.1</td><td>5</td><td>3.33</td><td>0–11.36</td><td>11.33</td><td>3</td></tr><tr><td>RAS5.2</td><td>5</td><td>112.88</td><td>92.95–137.95</td><td>45.01</td><td>307</td></tr><tr><td>RAS5.3</td><td>5</td><td>156.25</td><td>148.12–158.68</td><td>10.56</td><td>41</td></tr><tr><td>RAS6.1</td><td>6</td><td>110.41</td><td>102.08–118.41</td><td>16.33</td><td>16</td></tr><tr><td>RAS6.2</td><td>6</td><td>226.25</td><td>218.25–240.48</td><td>22.23</td><td>13</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Comparison of gene expression in crown root tips between ICML-IS 11139 and ICML-IS 11084. (A) Venn diagram representing the number of genes that are significantly differentially expressed according to three different statistical tests (EdgeR, DESeq and DESeq2; *p* value < 0.05) and the overlaps. (B) Main GO terms overrepresented in the DEG.</figcaption></figure>

Supplementary Table S2). Furthermore, the peak position of RAS5.1 on chromosome 5 is located 43 kb away from the SNP chr5_3282686 identified by GWAS (LFMM and EMMA). Likewise, the RAS5.3 spans through a genomic region containing two GWAS QTLs (LFMM, MLM and EMMA) located 113 kb and 244 kb away from the peak position of RAS5.3.

Altogether, the combination of GWAS and BSA analyses revealed genomic regions on chromosomes 5 and 6 controlling RAS/RT ratio in pearl millet.

**Comparison of gene expression in contrasted lines.** To further analyse the genes involved in rhizosphere formation, we compared gene expression in ICML-IS 11139 (low RAS/RT) and ICML-IS 11155 (high RAS/RT) roots. Production and secretion of root exudates occur along the root system³⁷, starting in the zone immediately behind the root tip³⁸. Similarly, root hair development occurs in the root tip. Thus, as these two processes seem to be the major determinants of root-soil aggregation in pearl millet, we hypothesized that genes controlling this trait might be preferentially expressed in the root tip. Phenotyping for RAS/RT ratio was performed at 28 days after planting, when the root system of pearl millet was made of one primary root and several crown roots possessing lateral roots³⁹. As crown roots make up most of the root system at this stage and to avoid noise due to sample heterogeneity (different root types), we therefore compared gene expression in the crown root tips (2 cm apex) of the two contrasted lines. RNAseq revealed 1270 genes with significant differences in gene expression between the two contrasted lines using three combined statistical tests (EdgeR, DESeq et DESeq2, *p* value < 0.05; Supporting Information Fig. 4). A gene ontology analysis on 742 genes with GO annotation out of the 1270 differentially expressed genes revealed a significant enrichment in GO terms associated with proteins
involved in molecular interactions (GO:0043531, ADP binding with lowest *p* value) and enzymatic reactions (GO:0016706, oxidoreductase activity for instance; Supporting Information Table S4).

We combined GWAS, BSA and gene expression analyses to identify candidate genes for RAS aggregation. We first focused our search for candidate genes in the QTL regions identified by GWAS that were coincident with regions of significance defined through BSA on chromosomes 5 and 6 (GWAS QTLs 5.1, 5.3, 5.5, 5.6 and 6.3). We assessed the annotated genes from the reference genome³ included in a 1 Mbp region centred around the most significant SNP position together with their expression data from the RNAseq experiment.

The most significant SNP marker in GWAS QTL 5.1 maps in chromosome 5 position 3,282,686 bp in an intergenic region between a cluster of four genes coding for glyoxylate reductase (*Pgl_GLEAN_10016760*, *Pgl_GLEAN_10016761*, *Pgl_GLEAN_10016762* and *Pgl_GLEAN_10016764*). Out of the four genes, one was differentially expressed in the contrasted lines for RAS/RT, the others showed a weak and variable expression level within the same genotype. Glyoxylate reductases are recycling enzymes that reduce glyoxylate to glycolate⁴⁰. Interestingly, the glyoxylate cycle plays an important role in the synthesis of malate, which is a major metabolite excreted in root exudates⁴¹. This region also contains *Pgl_GLEAN_10016765*, a gene coding for an arginase with significantly higher expression in ICML-IS 11139 (low RAS/RT). Arginases metabolise arginine and provide nitrogen for the synthesis of other essential amino acids during plant development and stress response mechanisms⁴². Large variations in arginine concentrations have been associated with changes in root exudate composition in plants exposed to drought⁴³.

The GWAS QTLs 5.5 and 5.6 are coincident with the same region of significance defined in BSA, RAS5.3. This 10.56 Mbp region contains 105 annotated genes in the reference genome. Interestingly, the most significant marker trait association for GWAS QTL 5.5 falls into a gene showing some homology to remorins (*Pgl_GLEAN_10037821*). Remorins are membrane proteins playing an important role in plant biotic interactions⁴⁴.

The most significant SNP in GWAS QTL 6.3 maps in chromosome 6 position 227,616,229 bp in a gene coding for a galactinol-sucrose galactosyltransferase that is more expressed in the low RAS/RT ratio line (*Pgl_GLEAN_10028942*). These are enzymes involved in the synthesis of raffinose⁴⁵, an oligosaccharide stored principally in seeds, roots and tubers. Accumulation of raffinose in wheat and tomato roots occurs in response to low P conditions⁴⁶,⁴⁷. Raffinose also accumulates in roots of pea seedlings exposed to water stress⁴⁸. In addition, the secretion of this oligosaccharide in root exudates is linked to the complex biotic interactions in the rhizosphere⁴⁹,⁵⁰.

We also looked for candidate genes associated with the most significant SNPs consistently identified by GWAS that do not coincide with regions of significance defined through BSA. We focused on GWAS QTLs accounting for large proportion of the phenotypic variance explained. The GWAS QTL 2.3 (PVE 12.4%) contains a cluster of five significant SNPs mapping in the same gene, *Pgl_GLEAN_10019483*, encoding an LRR receptor-like serine/threonine-protein kinase. This gene is strongly expressed in the root tip of both pearl millet lines. LRR receptor kinases are involved in the perception of signalling molecules⁵¹.

The GWAS QTL 6.2 (PVE 15.1%) consists of two SNPs markers on chromosome 6 mapping in an intergenic region between a cluster of four genes coding for acidic endochitinase (*Pgl_GLEAN_10020193*, *Pgl_GLEAN_10020194*, *Pgl_GLEAN_10020195* and *Pgl_GLEAN_10020196*), one of them with higher expression in the low aggregation line, the others with similar or weak expression in both lines. Endochitinase and chitinase-like proteins are defence related proteins with anti-fungal activity that are found in root exudates of different plant species⁵²⁻⁵⁴.

In chromosome 7, we found a group of five significant SNP markers within a 74 kb range defining the GWAS QTL 7.5 (PVE 14.4%). These SNPs were close to a gene encoding a putative chloroplastic dicarboxylate transporter that exchanges malate for succinate, fumarate and 2-oxoglutarate (*Pgl_GLEAN_10006630*). All these are important components of root exudates. The region also contains a gene coding for an ABC transporter G family member (*Pgl_GLEAN_10006636*). ABC transporters are involved in the transport of root exudates⁵⁵,⁵⁶.

## Discussion

Here, we investigated root system architectural traits with potential impact on rhizosheath formation in pearl millet. The presence of root hairs is essential for rhizosheath formation but the impact of root hair length and density on rhizosheath size varies considerably between plant species¹⁵. Studies in wheat and maize showed that root hair length is strongly correlated with rhizosheath weight²²,⁵⁷. Similarly, in foxtail millet, increased rhizosheath formation was found related with the plastic response in root hair formation (increases in root hair elongation and density) in dry soils¹⁷. This relationship was not as clear in crops such as barley¹⁸. However, a recent study shows an increased rhizosheath formation in barley grown in drying soil associated with auxin-promoted growth of root and root hairs as a consequence of ABA accumulation⁵⁸. In our study on pearl millet, we have identified a weak but significant correlation between rhizosheath formation, which is synonymous with root-adhering soil formation in this work, and root hair length (*p* = 0.005, *r*² = 0.077). This suggests that root hairs are involved in rhizosheath formation in pearl millet but that they play a limited role in our experimental conditions.

Root association with arbuscular mycorrhizal fungi (AMF) has been proposed to contribute to rhizosheath formation²⁰. In our study, we did not find any correlation between rhizosheath formation and AMF colonization rate in pearl millet suggesting that AMF colonization level is not an important driver in rhizosphere aggregation in this species.

Altogether, we hypothesise that rhizosheath formation or the aggregation of soil particles to the root in pearl millet is mainly driven by other traits. Root exudates and mucilaginous polymers released by root-associated microorganisms as well as the enzymatic activities linked to the crosstalk interactions occurring in the rhizosphere are prime candidates. Accordingly, the different orders of bacteria predominantly found in the rhizosphere of pearl millet lines with contrasted root soil aggregation suggests that the differences in rhizosheath formation
could be linked to crosstalks between the plant and microbial community³³,³⁴. Further work will be needed to test this hypothesis.

In the current study, the large variation in rhizosheath size in a genetically diverse group of inbred lines revealed a high heritability value for the trait ($H^2=0.72$). Although rhizosheath formation relies on a range of traits mainly related with root morphology and exudates, it has been found under genetic control in other cereal crops such as wheat²⁸,³¹ and barley¹⁸,⁵⁹, becoming a potentially interesting target trait for breeding¹³. Chromosome regions associated with rhizosheath size were identified in both crops, however few candidate genes underlying the QTL regions have been proposed. Interestingly, comparative evaluation of the multiple loci identified in these studies shows a lack of QTLs identified across diverse growing conditions suggesting, to some extent, a large QTL by environment interaction likely linked to the plasticity of rhizosheath formation.

Here, the combination of GWAS and BSA allowed the identification of five chromosome regions controlling rhizosheath size in pearl millet and ultimately some putative candidate genes based on gene annotations in the reference genome³. GWAS allowed the identification of 34 significant QTLs using the latent factor mixed model or LFMM³⁶ method. Many of these associations were confirmed using two other models for GWAS analysis (EMMA and MLM). The phenotypic variance explained (PVE) by these loci ranged from 11.2% to 14.7% suggesting that rhizosheath size as a complex trait determined by many QTLs of moderated effect in pearl millet. Consistently, studies in biparental and multiparental populations of wheat revealed several QTLs linked to rhizosheath formation with proportions of variation explained by QTLs around 5–10%²⁸,³¹. Nonetheless, one major QTL for the trait was also identified in wheat³¹.

Few genetic studies have identified genes potentially involved in rhizosheath formation and their predicted functions were mainly linked to root system morphogenesis and growth. For example, root hair length is a major driver determining rhizosheath size in wheat and, accordingly, genes coding for basic helix-loop-helix family of transcription factors that are known to control root hair development were identified as potential candidates underlying a rhizosheath QTL in that species²⁸. In barley, genes controlling cell division in root apical meristem at seedling stage and genes linked to tolerance to drought and cold were also identified as putative candidates underlying some genomic regions associated with rhizosheath size¹⁸.

In contrast, in our study, candidate genes were mostly related to plant metabolism and transport. Combining BSA and GWAS analyses revealed five co-localizing QTL regions. Candidate genes in these QTLs regions were mostly linked with root metabolic activities such as the synthesis of compounds commonly found in root exudates. For instance, the glyoxylate reductase and the arginase identified as putative candidates for QTL 5.1 are involved in the reduction and storage of essential compounds (i.e., glyoxylate and nitrogen) required for metabolic processes that mediate the synthesis of organic acids like malate and the synthesis of amino acids, respectively⁴²,⁶⁰. These are major primary metabolites of root exudates which variations in concentration can trigger plant adjustments to enhance root access and mobilisation of soil phosphate and nitrogen when these nutrients are limited⁶¹⁻⁶³. Further, these compounds have been found to promote chemotaxis of beneficial bacteria into the rhizosphere⁶⁴. In fact, a recent study showed how differences in malate concentration in root exudates impacted the composition of microbial communities associated with wheat and rice root systems⁶⁵. Another potential candidate gene identified for QTL 6.3, a galactinol-sucrose galactosyltransferase, is involved in the synthesis of raffinose, an oligosaccharide which variations in concentration in root exudates has been found to favour root colonisation by rhizosphere microbes⁵⁰.

Our genetic analysis is therefore fully consistent with our analysis showing that root architectural, root hair and AM symbiosis traits are not or poorly correlated with rhizosheath formation in pearl millet, and with our expression study that shows that genes involved in plant metabolism are differentially regulated between lines with contrasted rhizosheath size. It is also consistent with previous research showing differences in the rhizosphere metabolic activity of pearl millet lines with contrasted rhizosheath size³³. In this study, increased activity of enzymes such as chitinase and phosphomonoesterase was observed in the rhizosphere of pearl millet lines with larger rhizosheath (same contrasted RAS/RT lines used in the present work). We hypothesised that increased exudation in lines with larger rhizosheath size lead not only to an enhanced stability of root-adhering soil aggregates but also to a decrease of pH that could have stimulated these enzyme activities³³. Moreover, the amount of root exudate and the function of these enzymes could also impact the rhizosphere microbial communities promoting rhizosheath formation and explain the difference found in microbiota diversity in contrasted pearl millet lines³³. In conclusion, our physiological and genetic analysis suggest a central role for root exudation (quantitatively or qualitatively) in the regulation of rhizosheath formation in pearl millet. Rhizosheath formation seems to be controlled by many QTLs with small effects. We identified several candidate genes controlling this trait and future work will focus on the validation and characterization of the molecular mechanisms regulating rhizosheath formation in pearl millet. The rhizosheath has been proposed as a potential new target for breeding crops that are more tolerant to abiotic stresses¹³. Indeed, increased rhizosheath formation has been positively correlated with improved tolerance in drought stress in various cereals (wheat, barley, oat and foxtail millet)¹⁶⁻¹⁹ and in chickpea²⁷. Future work will also need to address the impact of rhizosheath formation in agronomical performances in pearl millet to better evaluate its potential for breeding.

## Materials and methods

**Plant materials.** A panel of 181 pearl millet inbred lines developed at the International Crops Research Institute for the Semi-Arid Tropics (ICRISAT, Niger) from landraces and improved open-pollinated cultivars representing the genetic diversity of the crop in West and Central Africa was used in this study³⁵.

Two inbred lines from this panel with contrasted rhizosheath size measured by the ratio between the mass of root-adhering soil (RAS) and root biomass (RT; RAS/RT ratio)³³: ICML-IS 11139 (small rhizosheath size
parent) and ICML-IS 11084 (large rhizosheath size parent) were selected for a bi-parental cross. The obtained F2 offspring was then used in a bulk segregant analysis (BSA).

Seeds were produced by ICRISAT at its Sadoré station (Niger) in compliance with relevant institutional, national, and international guidelines and legislation.

**Plant growth and measurement of soil aggregation.** Plants were grown for 28 days in “WM” shaped pots (WM 20-8-5, Thermoflan, Molières-Cavaillac) containing 1.5 kg of soil under natural light in a greenhouse in the ISRA/IRD Bel Air Campus in Dakar (Lat. 14.701778, Long. – 17.426229, altitude 9 m) as previously described³³.

For the GWAS analysis, pearl millet lines were sown according to a complete random block design with seven repetitions. Thinning was performed to have one plant per pot. Soil moisture was adjusted daily at water-holding capacity. Plant watering was stopped 24 h before harvesting to facilitate the separation of root-adhering soil (RAS) from bulk soil. Plants were harvested 28 days after planting by opening the pots gently and shaking the plant and its adhered soil at a constant speed (1100 rpm) for 1 min with a CAT S50 electric shaker (Cat Ingenieurbuero™) to separate the bulk soil from the RAS uniformly. Roots were then rinsed in a cup with demineralized water to collect RAS. The RAS was dried at 105 °C for three days and weighted. Roots and shoots were separated and dried at 65 °C for three days. The ratio between mass of RAS and mass of root tissue (RT; RAS/RT) was used to estimate the rhizosphere aggregation intensity (rhizosheath size) as previously described³⁴.

For BSA analysis, RAS/RT ratio was measured on 553 F2 individuals grown in five successive blocks of 119, 112, 131, 130 and 61 F2 plants. Each of these blocks included six replicates randomly distributed for each parental line. At the end of the experiment, leaf disk samples of 1.5 mm diameter were sampled for each individual plant and stored at – 80 °C for genotyping.

For correlation analyses between RAS/RT ratio and related root traits (root architecture, root hair length and density and interactions with AMF), 8 contrasting genotypes for rhizosheath size were analysed in 2018 and 2020 (n = 10 plants/genotype in 2018 and n = 6 plants/genotype in 2020).

**Root architecture.** Root architecture traits (length, average diameter, total root area) were measured using the WinRHIZO software (version 2012b) after scanning the roots using an Epson Perfection V700 scanner. Roots were separated in two groups based on their diameter according to Passot et al.³⁹: primary and crown roots (0.25 mm < diameters) and lateral roots (diameters < 0.25 mm).

Root hair length and density were measured on four plants per genotype using images of the root hair zone of three lateral roots per plant. Images were taken using an optical microscope (BX50F, Olympus) equipped with a digital camera (Micro Publisher 3.3 RTV). For each lateral root, the total number of root hairs was recorded over a distance of 0.5 mm using the Mesurim free software (http://acces.ens-lyon.fr/acces/logiciels/applications/mesurim/mesurim) and the length of 10 randomly selected root hairs was measured using the ImageJ software.

**Root colonization by arbuscular mycorrhizal fungi.** Intensity and frequency of root colonization by AMF were measured according to⁶⁶ after roots staining with Trypan blue following the method described by Phillips & Hayman⁶⁷. Stained root fragments were observed with a Nikon Labophot trinocular microscope. For each fragment, a score between zero and five was assigned according to the estimated proportion of root cortex colonized by AMF⁶⁶.

Frequency and intensity of root colonization were then computed using the following formulas:

$$ \text{Frequency (expressed in \%): } F = n/N \times 100 $$

where n is the number of fragments showing mycorrhizae and N, the number of observed fragments

$$ \text{Intensity (expressed in \%): } I = (95n5 + 70n4 + 30n3 + 5n2 + n1)/N $$

where n1, n2, n3, n4, n5 are the number of fragments scored respectively from 1 to 5 and N, the number of fragments observed.

**Heritability.** Broad sense heritability was computed using the following formula:

$$ H^2 = \frac{\text{Var(line)}}{\text{Var(line)} + \frac{\text{Var(res)}}{n_{\text{plant/line}}}} $$

where nplant/line is the average number of plants measured per line, Var(line) is the variance associated with lines and Var(res) is the residual variance.

Both variances are parameters of the following linear mixed model:

$$ \text{RSA} = \mu + \alpha_{\text{line}} + \varepsilon_{\text{res}}, $$

where μ is the overall mean soil aggregation, αline is the random effect attached to the lines with αline ~ N(0, Var(line)) and εres is the error term with εres ~ N(0, Var(res)).

**Genome wide association mapping and statistical analysis.** Genotyping by sequencing of this panel of inbred lines was previously reported³⁵. As a preliminary step, we used the genotypic matrix to estimate the population structure. Individual ancestry coefficients were estimated using the R package LEA v2.0⁶⁸. We
used a latent factor mixed model (LFMM) that considers ridge estimates and corrects for unobserved population cofounders, i.e. latent factors, to perform the GWAS³⁶. In addition, we ran the efficient mixed-model association (EMMA)⁶⁹ and mixed linear model (MLM) implemented in the R package GAPIT⁷⁰ to contrast the results. The proportion of the phenotypic variance explained by a QTL was determined by estimating the R² corrected for population structure of a linear model defined for the most significant SNPs.

**Bulk segregant analysis.** NGS-based BSA studies require establishing contrasted groups or bulks of lines to assess the differences in segregation of alleles using high-throughput sequencing⁷¹. The 10% extreme lines in the tails of the phenotype distribution for the RAS/RT ratio were selected and the corresponding leaf discs were pooled to form bulks of contrasted lines. Genomic DNA was isolated for each bulk using a MATAB (Mixed Alkyl Trimethyl Amonium Bromide) based method⁷² and enriched DNA libraries were constructed for which 32,860 predicted genes from the pearl millet reference genome³ were targeted using gene capture probes (myBaits®). High-throughput sequencing of the enriched DNA library was performed on an Illumina HiSeq platform by Novogene Company Limited (HK). Initial sequencing quality checks using FastQC version 0.11.5⁷³ were followed by trimming and quality filter steps on which adaptors, barcode sequences and low-quality reads (< 35 bp) were removed. Paired sequences were then retained and aligned to the pearl millet reference genome using the BWA MEM algorithm (BWA version 0.7.17-r1188)⁷⁴. Reads mapping at the target-enriched regions were used for SNP calling using the UnifiedGenotyper algorithm from GATK 3.7⁷⁵. Down-sampling limit (dcov) was increased from the default value of 250 to 9000 to ensure accounting for the maximum coverage reached at each position. Multi-allelic sites and those which exhibited a total allele frequency less than 0.25 were removed. In addition, sites with either low or high total sequencing depth (below the 25th and above the 95th percentiles respectively) were removed. SNPs with more than 50% missing data and minor allele frequency (MAF) under 5% were also excluded. Finally, the parental line ICML-IS 11139 (low RAS/RT ratio) was used as the reference genome for the cross to designate the alternate and reference SNP variants in the bulks.

Euclidean distance-based statistics⁷⁶ was used to measure the difference in allele frequency between the bulks. The Euclidean distance between allele frequencies of the bulks at each marker position (EDm) was calculated as follows:

$$ED_m = \sqrt{(f_{aL} - f_{AL})^2 + (f_{aH} - f_{AH})^2}$$

where f_a and f_A correspond to the allele frequency of the alternate and reference allele in the low bulk (L) and the high bulk (H) respectively.

In order to reduce the effect of sequencing noise and increase the signal of the differences in allele frequency, we then calculated the fourth power of the cumulative EDm value in windows of 100 consecutive markers⁷⁷,⁷⁸.

**Gene expression analyses.** Seeds from lines ICML-IS 11139 (low RAS/RT ratio) and ICML-IS 11155 (high RAS/RT ratio) were surface-sterilized and germinated in Petri dishes containing wet filter paper for 24 h in the dark at 27 °C. After two days, plants were transferred to hydroponic tanks containing liquid half Hoagland solution and grown for 15 days at 27 °C (12 h light/12 h dark). RNA was extracted from crown root tips (two cm apex) using the RNeasy Plant Mini Kit (QIAGEN). RNA-seq was performed by the MGX- Montpellier GenomiX Platform. Sequencing was performed on an Illumina HiSeq 2500. Three different statistical tests were used to identify differentially expressed genes: EdgeR⁷⁹, DESeq⁸⁰ and DESeq2⁸¹. GO terms enrichment was performed in the 1270 genes that were significantly differentially expressed between the two lines for all three statistical tests using the TopGO package in R.

**Statistical methods.** All statistical analyses were performed with R version 3.6.3⁸².

**Data availability**

The data that support the findings of this study are openly available at the National Center for Biotechnology Information (NCBI). Genotyping (GBS) data are available in genbank under reference number PRJNA492967 (GWAS) and PRJNA769524 (BSA). RNAseq data are available in the Gene Expression Omnibus (GEO) under reference GSE185425.

Received: 17 December 2021; Accepted: 13 May 2022

Published online: 02 June 2022

**References**

1. Burgarella, C. et al. A western Sahara centre of domestication inferred from pearl millet genomes. *Nat. Ecol. Evol.* 2, 1377–1380 (2018).

2. Debieu, M., Kanfany, G. & Laplaze, L. Pearl millet genome: Lessons from a tough crop. *Trends Plant Sci.* 22, 911–913 (2017).

3. Varshney, R. K. et al. Pearl millet genome sequence provides a resource to improve agronomic traits in arid environments. *Nat. Biotechnol.* 35, 969–976 (2017).

4. Lynch, J. P. Root phenotypes for improved nutrient capture: An underexploited opportunity for global agriculture. *New Phytol.* 223, 548–564 (2019).

5. York, L. M., Carminati, A., Mooney, S. J., Ritz, K. & Bennett, M. J. The holistic rhizosphere: Integrating zones, processes, and semantics in the soil influenced by roots. *J. Exp. Bot.* 67, 3629–3643 (2016).

6. Wissuwa, M., Mazzola, M. & Picard, C. Novel approaches in plant breeding for rhizosphere-related traits. *Plant Soil* 321, 409–430 (2009).
7. de la Fuente Cantó, C. et al. An extended root phenotype: the rhizosphere, its formation and impacts on plant fitness. *Plant J.* **103**, 951-964 (2020).

8. Turner, T., James, E. K. & Poole, P. The plant microbiome. *Genome Biol.* **14**, 209 (2013).

9. Chai, Y. N. & Schachtman, D. P. Root exudates impact plant performance under abiotic stress. *Trends Plant Sci.* https://doi.org/10.1016/j.tplants.2021.08.003 (2021).

10. Haichar, Z., Marol, C., Berge, O. & Rangel-castro, J. I. Plant host habitat and root exudates shape soil bacterial community structure. *ISME J.* **2**, 1221-1230 (2008).

11. Hinsinger, P., Bengough, A. G., Vetterlein, D. & Young, I. M. Rhizosphere: Biophysics, biogeochemistry and ecological relevance. *Plant Soil* **321**, 117-152 (2009).

12. Yang, J., Kloepper, J. W. & Ryu, C. M. Rhizosphere bacteria help plants tolerate abiotic stress. *Trends Plant Sci.* **14**, 1-4 (2009).

13. Ndour, P. M. S., Heulin, T., Achouak, W., Laplaze, L. & Cournac, L. The rhizosheath: From desert plants adaptation to crop breeding. *Plant Soil* **456**, 1-13 (2020).

14. Price, S. R. The roots of some North African dessert-glasses. *New Phytol.* **10**, 328-340 (1911).

15. Brown, L. K., George, T. S., Neugebauer, K. & White, P. J. The rhizosheath—A potential trait for future agricultural sustainability occurs in orders throughout the angiosperms. *Plant Soil* **418**, 115-128 (2017).

16. Basirat, M., Mousavi, S. M., Abbaszadeh, S., Ebrahimi, M. & Zarebanadkouki, M. The rhizosheath: A potential root trait helping plants to tolerate drought stress. *Plant Soil* **445**, 565-575 (2019).

17. Liu, T. Y. et al. Rhizosheath formation and involvement in foxtail millet (*Setaria italica*) root growth under drought stress. *J. Integr. Plant Biol.* **61**, 449-462 (2019).

18. George, T. S. et al. Understanding the genetic control and physiological traits associated with rhizosheath production by barley (*Hordeum vulgare*). *New Phytol.* **203**, 195-205 (2014).

19. Nambiar, E. Uptake of Zn65 from dry soil by plants. *Plant Soil* **44**, 267-271 (1976).

20. Pang, J., Ryan, M. H., Siddique, K. H. M. & Simpson, R. J. Unwrapping the rhizosheath. *Plant Soil* **418**, 129-139 (2017).

21. Moreno-Espindola, I. P., Rivera-Becerril, F., de Jesús Ferrara-Guerrero, M. & De León-González, F. Role of root-hairs and hyphae in adhesion of sand particles. *Soil Biol. Biochem.* **39**, 2520-2526 (2007).

22. Delhaize, E., James, R. A. & Ryan, P. R. Aluminium tolerance of root hairs underlies genotypic differences in rhizosheath size of wheat (*Triticum aestivum*) grown on acid soil. *New Phytol.* **195**, 609-619 (2012).

23. Galloway, A. F. et al. Cereal root exudates contain highly structurally complex polysaccharides with soil-binding properties. *Plant J.* **103**, 1666-1678 (2020).

24. Dakora, F. D. & Phillips, D. A. Root exudates as mediators of mineral acquisition in low-nutrient environments. *Plant Soil* **245**, 35-47 (2002).

25. Kolb, E., Legué, V., Auvergne, U. C. & Ferrand, C. Physical root-soil interactions. *Phys. Biol.* **14**, 065004 (2017).

26. Sasse, J., Martinoia, E. & Northen, T. Feed your friends: Do plant exudates shape the root microbiome?. *Trends Plant Sci.* **23**, 25-41 (2018).

27. Rabbi, S. M. F. et al. Plant roots redesign the rhizosphere to alter the three- dimensional physical architecture and water dynamics. *New Phytol.* **219**, 542-550 (2018).

28. Delhaize, E., Rathjen, T. M. & Cavanagh, C. R. The genetics of rhizosheath size in a multiparent mapping population of wheat. *J. Exp. Bot.* **66**, 4527-4536 (2015).

29. Mwafulirwa, L. et al. Barley genotype influences stabilization of rhizodeposition-derived C and soil organic matter mineralization. *Soil Biol. Biochem.* **95**, 60-69 (2016).

30. Mwafulirwa, L. et al. Genotypic variation in maize (*Zea mays*) influences rates of soil organic matter mineralization and gross nitrification. *New Phytol.* **231**, 2015-2028 (2021).

31. James, R. A. et al. Rhizosheaths on wheat grown in acid soils: Phosphorus acquisition efficiency and genetic control. *J. Exp. Bot.* **67**, 3709-3718 (2016).

32. Mwafulirwa, L. et al. Identification of barley genetic regions influencing plant-microbe interactions and carbon cycling in soil. *Plant Soil* **468**, 165-182 (2021).

33. Ndour, P. M. S. et al. Pearl millet genotype impacts microbial diversity and enzymatic activities in relation to root-adhering soil aggregation. *Plant Soil* **464**, 109-129 (2021).

34. Ndour, P. M. S. et al. Pearl millet genetic traits shape rhizobacterial diversity and modulate rhizosphere aggregation. *Front. Plant Sci.* **8**, 1288 (2017).

35. Debieu, M. et al. Response to early drought stress and identification of QTLs controlling biomass production under drought in pearl millet. *PLoS ONE* **13**, e0201635 (2018).

36. Caye, K., Jumentier, B., Lepeule, J. & François, O. LFMM 2: Fast and accurate inference of gene-environment associations in genome-wide studies. *Mol. Biol. Evol.* **36**, 852-860 (2019).

37. Haichar, Z., Santaella, C. & Heulin, T. Soil biology & biochemistry root exudates mediated interactions belowground. *Soil Biol. Biochem.* **77**, 69-80 (2014).

38. Schroth, M. N. & Snyder, W. C. Exudation patterns from bean seeds and hypocotyls and their effects on *Fusarium solanii* f. *phaseoli* (Abstr.). *Phytopathology* **52**, 751 (1962).

39. Passot, S. et al. Characterization of pearl millet root architecture and anatomy reveals three types of lateral roots. *Front. Plant Sci.* **7**, 1-11 (2016).

40. Hoover, G. J. et al. Characteristics of an *Arabidopsis* glyoxylate reductase: general biochemical properties and substrate specificity for the recombinant protein, and developmental expression and implications for glyoxylate and succinic semialdehyde metabolism in planta. *Botany* **85**, 883-895 (2007).

41. Fernie, A. R. & Martinoia, E. Malate. Jack of all trades or master of a few?. *Phytochemistry* **70**, 828-832 (2009).

42. Siddappa, S. & Marathe, G. K. What we know about plant arginases?. *Plant Physiol. Biochem.* **156**, 600-610 (2020).

43. Gargallo-Garriga, A. et al. Root exudate metabolomes change under drought and show limited capacity for recovery. *Sci. Rep.* **8**, 1-15 (2018).

44. Jarsch, I. K. & Ott, T. Perspectives on remorin proteins, membrane rafts, and their role during plant-microbe interactions. *Mol. Plant-Microbe Interact.* **24**, 7-12 (2011).

45. Lehle, L. & Tanner, W. The function of myo-inositol in the biosynthesis of raffinose. *Eur. J. Biochem.* **38**, 103-110 (1973).

46. Nguyen, V. L., Palmer, L., Roessner, U. & Stangoulis, J. Genotypic variation in the root and shoot metabolite profiles of wheat (*Triticum aestivum* L.) indicate sustained, preferential carbon allocation as a potential mechanism in phosphorus efficiency. *Front. Plant Sci.* **10**, 1-14 (2019).

47. Sung, J. et al. Plant Science Metabolomic profiling from leaves and roots of tomato (*Solanum lycopersicum* L.) plants grown under nitrogen, phosphorus or potassium-deficient condition. *Plant Sci.* **241**, 55-64 (2015).

48. Lahuta, L. B., Pluskota, W. E., Stelmaszewska, J. & Szabli, J. Dehydration induces expression of GALACTINOL SYNTHASE and RAFFINOSE SYNTHASE in seedlings of pea (*Pisum sativum* L.). *J. Plant Physiol.* **171**, 1306-1314 (2014).

49. Fang, W. & Leger, R. J. S. Mrt, a gene unique to fungi, encodes an oligosaccharide transporter and facilitates rhizosphere competency in *Metarhizium robertsii*. *Plant Physiol.* **154**, 1549-1557 (2010).

50. Liu, Y. et al. Identification of root-secreted compounds involved in the communication between cucumber, the beneficial *Bacillus amyloliquefaciens*, and the soil-borne pathogen *Fusarium oxysporum*. *Mol. Plant-Microbe Interact.* **30**, 53-62 (2017).
51. Chakraborty, S., Nguyen, B., Wasti, S. D. & Xu, G. Plant leucine-rich repeat receptor kinase (LRR-RK): Structure, ligand perception, and activation mechanism. *Molecules* **24**, 3081 (2019).

52. De-la-Peña, C. et al. Root secretion of defense-related proteins is development-dependent and correlated with flowering time. *J. Biol. Chem.* **285**, 30654–30665 (2010).

53. Nóbrega, F. M., Santos, I. S., Da Cunha, M., Carvalho, A. O. & Gomes, V. M. Antimicrobial proteins from cowpea root exudates: Inhibitory activity against *Fusarium oxysporum* and purification of a chitinase-like protein. *Plant Soil* **272**, 223–232 (2005).

54. Tesfaye, M., Denton, M. D., Samac, D. A. & Vance, C. P. Transgenic alfalfa secretes a fungal endochitinase protein to the rhizosphere. *Plant Soil* **269**, 233–243 (2005).

55. Badri, D. V. et al. An ABC transporter mutation alters root exudation of phytochemicals that provoke an overhaul of natural soil microbiota. *Plant Physiol.* **151**, 2006–2017 (2009).

56. Baetz, U. & Martinoia, E. Root exudates: The hidden part of plant defense. *Trends Plant Sci.* **19**, 90–98 (2014).

57. Adu, M. O. et al. Quantifying variations in rhizosheath and root system phenotypes of landraces and improved varieties of juvenile maize. *Rhizosphere* **3**, 29–39 (2017).

58. Zhang, Y. et al. Abscisic acid mediates barley rhizosheath formation under mild soil drying by promoting root hair growth and auxin response. *Plant. Cell Environ.* **44**, 1935–1945 (2021).

59. Gong, X. & McDonald, G. QTL mapping of root traits in phosphorus-deficient soils reveals important genomic regions for improving NDVI and grain yield in barley. *Theor. Appl. Genet.* **130**, 1885–1902 (2017).

60. Igamberdiev, A. U. & Eprintsev, A. T. Organic acids: The pools of fixed carbon involved in redox regulation and energy balance in higher plants. *Front. Plant Sci.* **7**, 1042 (2016).

61. Carvalhais, L. C. et al. Root exudation of sugars, amino acids, and organic acids by maize as affected by nitrogen, phosphorus, potassium, and iron deficiency. *J. Plant Nutr. Soil Sci.* **174**, 3–11 (2011).

62. Mora-Macias, J. et al. Malate-dependent Fe accumulation is a critical checkpoint in the root developmental response to low phosphate. *Proc. Natl. Acad. Sci. USA* **114**, E3563–E3572 (2017).

63. Canarini, A., Kaiser, C., Merchant, A., Richter, A. & Wánek, W. Root exudation of primary metabolites: Mechanisms and their roles in plant responses to environmental stimuli. *Front. Plant Sci.* **10**, 157 (2019).

64. Feng, H. et al. Identification of chemotaxis compounds in root exudates and their sensing chemoreceptors in plant-growth-promoting rhizobacteria *Bacillus amyloliquefaciens* SQR9. *Mol. Plant-Microbe Interact.* **31**, 995–1005 (2018).

65. Kawasaki, A. et al. Manipulating exudate composition from root apices shapes the microbiome throughout the root system. *Plant Physiol.* **1**–17 (2021).

66. Trouvelot, A., Kough, J. L. & Gianinazzi-Pearson, V. Mesure du taux de mycorhization VA d'un système radiculaire. Recherche de méthode d'estimation ayant une signification fonctionnelle. in *The Mycorrhizae, Physiology and Genetics* (eds. Gianinazzi-Pearson, V. & Gianinazzi, S.) 217–221 (INRA Presse, 1986).

67. Phillips, J. M. & Hayman, D. S. Improved procedures for clearing roots and staining parasitic and vesicular-arbuscular mycorrhizal fungi for rapid assessment of infection. *Trans. Br. Mycol. Soc.* **55**, 158–161 (1970).

68. Frichot, E. & François, O. LEA: An R package for landscape and ecological association studies. *Methods Ecol. Evol.* **6**, 925–929 (2015).

69. Kang, H. M. et al. Efficient Control of Population Structure in Model Organism Association Mapping. *1723*, 1709–1723 (2008).

70. Lipka, A. E. et al. GAPIT: Genome association and prediction integrated tool. *28*, 2397–2399 (2012).

71. Nguyen, K. L., Grondin, A., Courtois, B. & Gantet, P. Next-generation sequencing accelerates crop gene discovery. *Trends Plant Sci.* **24**, 263–274 (2019).

72. Mariac, C. et al. Diversity of wild and cultivated pearl millet accessions (*Pennisetum glaucum* [L.] R. Br.) in Niger assessed by microsatellite markers. *Theor. Appl. Genet.* **114**, 49–58 (2006).

73. Andrews, S. FastQC: a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc (2010).

74. Li, H. & Durbin, R. Fast and accurate short read alignment with burrows-wheeler transform. *Bioinformatics* **25**, 1754–1760 (2009).

75. McKenna, A. et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. *Genome Res.* **20**, 1297–1303 (2010).

76. Hill, J. T. et al. MMAPPR: Mutation mapping analysis pipeline for pooled RNA-seq. *Genome Res.* **23**, 687–697 (2013).

77. Omboki, R. B. et al. Pooled mapping of quantitative trait loci conferring male sterility—Conditioned glume split in rice (*Oryza sativa*). *Plant Breed.* **137**, 848–856 (2018).

78. Zhang, H. et al. QTG-Seq accelerates QTL fine mapping through QTL partitioning and whole-genome sequencing of bulked segregant samples. *Mol. Plant* **12**, 426–437 (2019).

79. Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR : a Bioconductor package for differential expression analysis of digital gene expression data. **26**, 139–140 (2010).

80. Anders, S. & Huber, W. Differential expression analysis for sequence count data. *Genome Biol.* **11**, R106 (2010).

81. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biol.* **15**, (2014).

82. R. core team (2018). R: A Language and Environment for Statistical Computing. https://www.r-project.org.

## Acknowledgements

This work was supported by the French National Research Institute for Sustainable Development (IRD), the French Agence Nationale de la Recherche (ANR grant RootAdapt ANR17-CE20-0022-01 to LL), the CGIAR Research Programme on Grain Legumes and Dryland Cereals (GLDC), the NewPearl grant in the frame of the CERES initiative by the Agropolis Fondation (AF 1301-015 to LL as part of the ‘Investissement d’avenir’ ANR- I0-LABX-0001-01 under the frame of I-SITE MUSE ANR-16-IDEX-0006), and by the Fondazione Cariplo (No FC 2013-0891).

## Author contributions

Y.V., L.C. and L.L. conceptualized and supervised the research; L.L. acquired the funding. C.F.C., M.N.D., P.M.S.N., M.D., A.G., S.P., C.B., and M.P. carried out the measurements and performed formal analyses. All authors discussed and evaluated the data. C.F.C., M.N.D., P.M.S.N., A.G., Y.V., L.C. and L.L. wrote the first draft of the manuscript; all authors revised the manuscript and gave final approval for publication.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-022-13234-w.

**Correspondence** and requests for materials should be addressed to L.L.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_11.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022