Set up from the beginning: The origin and early development of cassava storage roots

Abstract Despite the importance of storage root (SR) organs for cassava and the other root crops yield, their developmental origin is poorly understood. Here we use multiple approaches to shed light on the initial stages of root development demonstrating that SR and fibrous roots (FR) follow different rhizogenic processes. Transcriptome analysis carried out on roots collected before, during and after root bulking highlighted early and specific activation of a number of functions essential for root swelling and identified root‐specific genes able to effectively discriminate emerging FR and SR. Starch and sugars start to accumulate at a higher rate in SR before they swell but only after parenchyma tissue has been produced. Finally, using non‐destructive computed tomography measurements, we show that SR (but not FR) contain, since their emergence from the stem, an inner channel structure in continuity with the stem secondary xylem, indicating that SR derive from a distinct rhizogenic process compared with FR.


| INTRODUCTION
Producing enough food, in a sustainable manner, to meet the needs of an increasing global population is one of the greatest challenges we are confronted with nowadays. Cassava (Manihot esculenta) is one of the most important staple food crops providing nutritional security to over 800 million people worldwide. With its storage roots (SR) consisting mostly of starch (up to 90% dry weight; Montagnac et al., 2009), and its resilience to heat and drought stresses, cassava can efficiently respond to critical present and future climate change challenges. Despite 6000 years of cassava domestication (Bredeson et al., 2016) and the relevance of its roots as staple food for millions of people in tropical and subtropical areas, relatively little is currently known about the mechanism of SR formation and starch accumulation of this important crop, and of root crops in general (Hershey et al., 2012;Jammer et al., 2020;Villordon et al., 2014).
Tuberous root development is a complex process involving intense crosstalk among phytohormones, metabolites and regulation of gene transcription to manage cell division and elongation for radial growth and starch accumulation (Li et al., 2016;Utsumi et al., 2020).
The anatomical process of SR formation has been well described in classical studies at the light microscopy resolution level (Kokubu, 1973;McCormick, 1916;Togari, 1950;Wilson & Lowe, 1973). Higher activity of the vascular cambium, achieved through cell divisions and expansions in the primary and secondary meristems, lead to an expansion of the root diameter (Togari, 1950), which exhibits a lower degree of lignification in the stele. First, there is rapid division of cambium cells along the root axis. This is followed by vascular cambium production of numerous parenchymatous tissues, resulting in the thickening of the root diameter needed to accumulate large quantities of starch (Naconsie et al., 2016).
Based on the observation that cassava SR start to swell in a bunch of fine-branched apparently undifferentiated fibrous roots (FR), they are assumed to develop from FR upon receiving one or more specific, but unknown signals (Alves, 2009;El-Sharkawy, 2004).
Many questions remain open regarding this hypothesis, for example, how and why certain FR are selected to become SR, and which (and when) the signal(s) are delivered to (or activated in) the bulking roots.
The mechanism of storage-organ growth has recently been investigated in various root crops using different molecular, microscopy and metabolite measurement technologies. Upregulation of sucrose (Zhang et al., 2017), starch and storage protein biosynthesis, accompanied by delignification (Firon et al., 2013;Sun et al., 2015), cell wall and tissue functional differentiation (Ding et al., 2020;Firon et al., 2013;Kuznetsova et al., 2020) and phytohormone activation (Sojikul et al., 2015) are the major modifications identified in the thickening transition stage and final development into SR. Thus, even though differences between FR and SR have been identified, the earliest stages of SR development have not been studied, particularly at the molecular level, in any of the root crops, and the hypothetical signals changing an FR into an SR have yet to be discovered.
In this study, we used systematic approaches to study the functional differentiation between FR and SR. We provide evidence of different genetic, morphologic and metabolic behaviours of FR and SR starting from their emergence from a planting stem and, therefore, earlier before SR swelling. Non-destructive X-ray computed tomography (CT) analysis confirmed the morphological differences during the root emergence stages. Understanding the regulatory processes and the underlying genetic components controlling root growth, and thereby root system architecture, is not only a fascinating question of basic plant biology but also essential for understanding cassava and all other root crops productivity, and is key to breed and engineer better-performing plants with a potential to feed the future world population in the face of global climate change.

| Planting and sample materials
The study was conducted at the International Institute of Tropical Agriculture in Ibadan, Nigeria. Freshly M. esculenta (cassava) cut stems of the TMEB419 genotype as well as KALESO, TMEB14, TMS-IBA090576 genotypes, obtained from the IITA Genebank collection, were used as planting material for the study. Trials were laid out in a randomized complete block design with three replications. Healthy stem cuttings, each 25 cm in length were vertically planted in a flat seedbed at a spacing of 1 × 1 m. Samples of FR, potential FR (PFR), SR, potential SR (PSR), source leaves (SOL) and sink leaves (SIL) were obtained from single or multiple organs as further specified. Samples were separated in different biological replicates, cut in small pieces, mixed and snap-frozen in liquid N 2 .
2.2 | RNA extraction and quantitative real-time PCR (RT-qPCR) RNA was extracted from samples of SIL, SOL and both PFR/FR and PSR/SR tissue by combining cetyl trimethylammonium bromide (CTAB)-extraction method and spin-column-based purification methods. Total nucleic acid was extracted by using the CTAB (2% CTAB, 2% PVP-40, 20 mM Tris-HCl, pH 8.0, 1.4 M NaCl, 20 mM ethylenediaminetetraacetic acid) extraction protocol as described by Li et al. (2008) with some modifications. Approximately 500 mg samples were ground in liquid nitrogen and mixed with 1 ml pre-heated CTAB extraction buffer. Samples were incubated at 65°C for 15 min with intermittent vortexing, then subjected to centrifugation (16 000 g at 4°C for 5 min). The supernatant was mixed with an equal volume of cold chloroform: isoamyl alcohol (24:1) before centrifugation (16 000 g at 4°C for 10 min). The supernatant was added to 0.6 volumes of cold isopropanol and gently mixed to precipitate the nucleic acids. The pellet was washed with 70% ethanol, air-dried and dissolved in nuclease-free water. After DNaseI treatment, the resulting RNA was cleaned up using the kit RNA Clean & Concentrator™ (Zymo Research) according to the manufacturer's instructions.
The RNA concentration in the samples was normalized based on Qubit assays (Invitrogen, Thermo Scientific) and gel observations, and 1 µg of quality-confirmed RNA was subjected to complementary DNA (cDNA) synthesis as follows. Reverse transcription was completed with 1 µg of total RNA in 20 µl. First-strand cDNA synthesis was performed using the M-MuLV Reverse Transcriptase (New England Biolabs) according to the manufacturer's instructions.
The reaction contained 1 µg of RNA, 2 µl of 10× RT random primer and 1 µl of 10 mM dNTPs. The mixture was denatured at 65°C for 5 min and placed on ice immediately. Next, 4.4 µl of the RT master mix (consisting of 2 µl 10× RT Buffer, 1 µl of M-MuLV RT, 0.2 µl of RNasin) were added and the mixture incubated at 25°C for 5 min, at 42°C for 1 h, and 65°C for 20 min. The cDNA synthesized served as a template for qPCR amplification of specific target genes.
Briefly, qPCR reactions containing 10 ng of cDNA, 400 nM of each forward and reverse primers and Luna ® Universal qPCR Master Mix (New England Biolabs ® ) were combined in a total volume of 12.5 µl. Three technical replicates per sample were used except for single plants (Figure 7) where we used four technical replicates. The reactions were conducted in a LightCycler ® 480 Instrument (Roche) using the following cycling profile: 10 s at 95°C, followed by 40 cycles of 3 s at 95°C and 30 s at 60°C. Data obtained were converted into relative gene expression using the 2 C -ΔΔ t method corrected for the starch granules were coated with platinum palladium and imaged with a Hitachi SU5000 SEM using a 5 kV electron beam. 3 | RESULTS

| SR development
Cassava roots, vegetatively propagated from stem cuttings (stakes) in field production systems, develop with a certain variability in time and consistency from stake to stake ( Figure S1). Nevertheless, depending on the genotype and the environmental conditions, SR, arising from both nodes and the basal cut, become evident starting from 6 to 8 weeks after planting (wap, Figure 1a). At early developmental stages, the emerging and elongating roots appear morphologically similar. However, in some cases and upon careful observation, we noticed two different root morphologies never described before: (i) thinner, softer and lighter in colour and (ii) thicker, stiffer and darker (Figures 1b and S1). The morphological difference becomes more visible over time (between 4 and 8 wap, Figure 1d), albeit it is faintly discernible even shortly after root emergence from the stake ( Figure S1). We named the first type 'potential fibrous roots' (PFR) and the second type 'potential storage roots' (PSR) and investigated whether the fate of PSR and PFR was indeed to develop into SR and FR, respectively. To this aim, we Putatively morphologically different PSR and PFR (TP1 and TP3) and FR and SR (TP7) were collected, and the different tissue samples were used to investigate gene expression and metabolic changes occurring among them over time.

| Functional gene expression analysis
To To identify functions and processes specifically related to SR formation, we selected genes significantly (p < 0.05) differentially expressed between PFR/FR and PSR/SR for each of the three TP and used GO to cluster them (Tulipano et al., 2007). The obtained clusters contain genes with similar functionalities identified by similar GO terms and are ranked based on changes of GO representation between time points in each GO category ( Figure 2 and Table S1).
The expression fold change of the genes associated to each functional cluster was compared over the three TPs (Figure 2, heatmaps).  Table S1). The expression of most of the genes relevant to the functions regulated in PSR between TP1 and TP3, remains largely unchanged until TP7, after root bulking (Figure 2, heatmaps).
On the other hand, the expression of another set of genes controlling functions and processes implicated in cell growth and differentiation (i.e., protein folding, cytoprotective unfolded protein response, amino acid synthesis and intercellular exchanges) and specifically regulated at TP3 in PSR becomes inversely regulated in SR at TP7, when root swell is appreciable ( Figure 2 and Table S1).
Root thickening, mainly driven by increased production of vascular and parenchymatous tissues, relies on the activation of a number of specific genes. In particular, vascular formation is strongly regulated by the transcription factor gene family VASCULAR-RELATED NAC-DOMAIN PROTEIN 1-7 (VND1-7). We observed that  (Kubo, 2005;Zhao et al., 2007), are both up-regulated in cassava PSR/SR relative to PFR/FR (Table S2). On the other hand, genes regulating ectopic secondary cell wall formation in various tissues, such as NAC SECONDARY WALL THICKENING PROMOTING FACTORs (NST1 and NST2) (Ko et al., 2007;Mitsuda et al., 2005;Zhong et al., 2006),  (Table S2).
As the main function of SR is starch accumulation, we monitored the regulation of the major genes supporting starch metabolism (Pfister & Zeeman, 2016)  probably not stored in the vacuole in SR, but rather directed towards starch biosynthesis or provide energy for root growth (Figure 3 and Table S3).
Downstream of SWEET, genes encoding sucrose synthases (SuSy, which catabolise imported sucrose), Glucose 6-phosphate (Glc6P)/phosphate translocators (GPTs, which import glucose 6-phosphate into amyloplasts), and ADP glucose pyrophosphorylase (AGPase, which synthesizes the substrate for starch synthesis) were up-regulated in PSR compared with PFR at TP3 (and in most cases at TP7; Figure 3). We did not detect up-regulation of cytosolic invertase (cytINV, which represents an alternative route for sucrose catabolism). Importantly, genes encoding Soluble Starch Synthase (SSSs) and Branching Enzymes (BEs)-both crucial for starch synthesis-were also up-regulated in PSR/SR compared with PFR/FR at the later time points (Figure 3 and Table S3). Key genes involved in starch degradation such as β-amylases (BAMs; which hydrolyse starch polymers to maltose) and pGlcT (which exports glucose derived from starch to the cytosol) (Lloyd & Kossmann, 2015) were down-regulated in PSR/SR compared with PFR/FR. Interestingly, however, disproportionating enzyme (DPE1, which metabolises malto-oligosaccharides) was up-regulated. This enzyme is involved in starch degradation but may also play a role in starch biosynthesis since malto-oligosaccharides are required to initiate starch polymer biosynthesis, and are also by-products of the biosynthesis process (Ball & Morell, 2003;Colleoni et al., 1999;Tetlow & Emes, 2014).
This overtime gene expression comparison between PFR/FR and PSR/SR indicates that a clear activation of starch production in PSR starts at TP3 while genes involved in the differentiation of PSR/SRspecific tissues are already active at TP1.

| Root-specific gene expression
To gain more insight into the transcriptomic signature of PSR/SR and PFR/FR and to identify potential root-specific genes, we extracted, from the total gene expression data, genes with FPKM ≤ 3 in both leaf tissues (SIL and SOL) at all time points and considered them as genes not expressed in leaves. From this first selection, genes with an FPKM ≥ 5 in at least one root tissue (PFR/FR or PSR/SR) and in at least one time point (TP1, TP3, or TP7) were picked out as candidate root-specific genes to be analysed for differential expression between PSR/SR and PFR/FR.
Based on these conditions, we identified five genes, all involved in cell growth and development, that were highly expressed in PFR/ FR and lower in PSR/SR (Figure 4a and Table S4

| FR and SR expression marker genes
We sought to confirm the differential expression behaviour of the root-specific genes identified and to investigate their potential use as expression markers specific for PFR/FR or PSR/SR in the initial stages of development. Therefore, we selected seven of the root-  Figure 5).
We tested whether these genes show similar expression behaviour also in other cassava genotypes, which is crucial to use them as molecular markers and potential biotechnological targets for crop improvement. The differential expression of FR-UP1, SR-UP3,  Figure S2). These three genotypes feature high SR dry matter content (%dry weight/fresh weight) but contrasting performances in root size (RS), root number (RN) and root weight (RW) traits ( Figure S2).  (Figure 6a).

| Metabolites analysis
At TP7, when SR are visibly swollen, and FR continue elongating, the soluble sugars content remained stable or slightly decreased in FR (negative R 2 ). In contrast, fructose, glucose and sucrose content increased in SR relative to TP1 and to FR (Figure 6a). This increase in soluble sugars presumably provides energy for continued PSR expansion and for starch biosynthesis in the bulking SR.
As expected, the magnitude of starch concentration correlates with SR swelling (Figure 6a). Starch is present in very low concentration in both PFR and PSR at TP1 and remained low in PFR/FR over time. In PSR/SR, starch content started to increase at TP3 and in bulked roots (TP7), starch had increased 40-fold relative to FR (Figure 6a).
We used scanning electron microscopy (SEM) to visualise starch granule morphology in PFR/FR and PSR/SR over the different stages of root development. While granule shape appeared similar among all samples (Figure 6b), the granules size increased over time, with massive increases in PSR/SR (Figure 6c). In PFR/FR, most granules at all three TPs were small (~2 µm), but the frequency of bigger granules increased at the latter TPs (Figure 6d). In PSR/SR, the granules diameter significantly increases over time (Figure 6e). At TP1, all granules were small, as in PFR/FR. At TP3 a significant fraction of the granules had diameters >3 µm. Interestingly, at TP7, two populations of granules could be identified: smaller granules with diameters between 1 and 5 µm and granules with diameters between 5 and 10 µm (Figure 6e).

| Gene expression and carbohydrate levels in single roots of single plants
Our results suggest that sugar levels do not change significantly between PFR and PSR before root bulking (TP7) whereas differential gene expression in the two organs occurs already at TP1, far earlier than the onset of morphological organ changes.
Seeking further confirmation that PSR and PFR differ early in their development, we selected two genes (FR-UP1 and SR-UP3) having contrasting levels of expression in PSR and PFR ( Figure 5). We

| Non-destructive analysis of cassava root apparatus development
To gain more insight into the development process of the cassava root apparatus, non-destructive X-ray CT measurements were taken over a 13-week time course in pot-grown plants. Figure Figure S4), while in all swell SR the channel structure was detectable since emergence and enlarged continuously with root thickening (Figures 9 and S3). However, not all roots that displayed an inner channel swelled into an SR, which is consistent with our observation that not all roots visible scored as PSR swelled, despite being thicker and darker than roots scored as PFR. Interestingly, upon cutting of all formed SR from the planted stem, swelling reinitiated in the roots containing channel structure (PSR) at a belowground growth rate similar to that measured before SR cutting ( Figure S5). When both SR and PSR (i.e., all channel-containing roots) were physically cut off, we  Figure S3. Green arrows point at the low-density channel. In the images, light colour corresponds to high density. Scale bars, 10 mm. FR, fibrous roots; PFR, potential fibrous roots; PSR, potential storage roots; SR, storage roots [Color figure can be viewed at wileyonlinelibrary.com] measured significantly longer time until plant roots started bulking again, and it occurred at an initial growth rate lower than before the cut ( Figure S5). We never observed SR developing from PFR/ FR left attached to the stem; therefore, bulking initiated from newly emerged PSR. These results indicate that PSR have the potential to start cell proliferation and starch accumulation so as to turn into SR, and in their absence PFR/FR (roots without channel structure) will not turn into SR. some of them undergo the process of starch accumulation probably triggered by a signal yet to be uncovered.

| DISCUSSION
Our observation of contiguous connection of the root and stem secondary xylem, previously only reported in swollen SR (Chaweewan & Taylor, 2015), is consistent with a distinct rhizogenesis of PSR, possibly involving a direct extension of the secondary xylem tissue of the stem, which itself accumulates starch (up to 42% of dry mass, Wei et al., 2015). Therefore, the presence of stemderived secondary xylem in PSR and not in PFR could explain why the former root type is preset to immediately activate tissue differentiation and ultimately starch accumulation.
Our results also shed light on the debate over where on the planting stake SR origin. Our and previous observations (Lowe & Mahon, 1982) of SR arising from both nodes and (prevalently) the basal cut of stems planted either in the field or in large pots (≧30 cm diameter) definitely rebut the hypothesis that cassava SR only form at stem nodes (Chaweewan & Taylor, 2015). However, yield reduction observed upon increasing cassava planting density (Silva et al., 2013) indicates that space limitation hinders SR development, and suggests that the sideways growth of PSR emerging from nodes is favoured in small pots (12 cm) as compared with basal roots (Chaweewan & Taylor, 2015).
We illustrated that PSR of field cultivated plants can be occasionally discriminated from PFR by a darker and stiffer appearance becoming progressively visible over time, conferred by the accumulation of suberin, phenolic compounds and other metabolites when periderm replaces epidermis during root secondary growth (Campilho et al., 2020). This process must be activated earlier and more intensively in PSR compared with PFR. The downregulation of genes associated with lignification and up-regulation of genes involved in differentiation of metaxylem-like vessels, detected in young PSR, suggest that the thin channel structure measured by CT in PSR (and not in PFR) corresponds to the early establishment of xylem parenchyma cells, which only later become distinguishable as hallmarks of SR via staining and histological analysis (Chaweewan & Taylor, 2015).
An important new insight provided by our study is that differential expression of selected genes can effectively discriminate between PFR and PSR in several cassava genotypes already 1 wap. The availability of early selection markers will enable correct discrimination and comparison between the two root types and will be critical to shed additional light on the process of SR development and to accelerate breeding selection. While orthologues of these marker genes are generally involved in growth and development, except for Pt2L4 (SR-UP4), P54 (SR-UP5) and Knotted 1-like (SR-UP2), the role of the selected genes in cassava SR development and their potential application in biotechnological approaches to engineer high yielding