Auxin signaling and vascular cambium formation enable storage metabolism in cassava tuberous roots

Auxin-mediated activation of secondary growth and subsequent KNOX/BEL expression coincide with active storage metabolism in xylem parenchyma cells of the cassava tuberous root.


Introduction
Cassava (Manihot esculenta) is an agronomically important root crop species grown in tropical and subtropical regions of the world (Carvalho et al., 2018). Over 60% of the global annual cassava yield is produced in sub-Saharan Africa (http:// faostat3.fao.org/), despite the crop being almost exclusively grown by smallholder farmers with limited resources. While grown on a large total field area, the yield per hectare is often low due to suboptimal agronomic practice, the lack of fertilizer, and/or high pathogen incidence. In addition, progress in cassava breeding lags behind that of other crops such as maize, wheat, or rice. Current cassava breeding focuses mostly on resistance and complex yield traits, such as dry matter content and root yield. A better understanding of cassava storage root development could facilitate breeding or biotechnological efforts aimed at improved storage root yield.
Cassava is clonally propagated through the planting of stem cuttings from which new shoots and roots can emerge. Cassava storage roots develop from stem-derived roots through the formation of a vascular cambium and subsequent secondary root growth (Chaweewan and Taylor, 2015;Mehdi et al., 2019). The molecular mechanisms of secondary growth have mainly been studied in the model systems Arabidopsis thaliana and poplar (Populus trichocarpa), and, while many details are still elusive, much has been learned in the recent past (Fischer et al., 2019). The vascular cambium develops from a procambium in between primary phloem and xylem tissue, and controls the formation of secondary phloem and xylem (Ragni and Greb, 2018;Fischer et al., 2019). The initiation and maintenance of meristematic stem cells is largely regulated by the phytohormone auxin (Snow, 1935;Björklund et al., 2007;Agusti et al., 2011). Auxin is synthesized in the shoot apex and moves into the root via a polar transport mechanism governed by auxin influx (AUX) and efflux carriers (PIN) (Bennett et al., 1996;Carraro et al., 2012;Robert et al., 2015). Here it leads to the formation of a stem cell organizer via CLASS III HOMEODOMAIN LEUCINE ZIPPER (HD-ZIPIII) transcription factor activity that regulates the position and maintenance of the vascular cambium (Smetana et al., 2019). Auxin is not the only phytohormone involved in vascular cambium development. Cytokinin and gibberellin (GA) have a complex interplay with auxin in the meristem of poplar stems (Immanen et al., 2016). Cytokinin plays a major role in defining cell identities and controlling cell division, while GA promotes cell elongation and the deposition of secondary cell walls and lignification (Eriksson et al., 2000;Biemelt et al., 2004;Fischer et al., 2019).
While some studies have recently indicated the importance of early auxin signaling events, as well as reduced GA signaling during storage root formation in sweet potato (Noh et al., 2010;Firon et al., 2013;Singh et al., 2019), less is known about storage root initiation in cassava. The importance of auxin and GA, as well as other phytohormones, could be shown for in vitro root formation (Sojikul et al., 2015;Utsumi et al., 2020). In addition, Ding et al. (2020) recently analyzed the impact of transcriptional and post-transcriptional regulation of starch and cell wall metabolism as well as hormone response genes in already established cassava storage roots. However, processes during the early stages of storage root development remain to be elucidated in cassava.
In this study, we focused on the morphological and transcriptional changes occurring during the early stages of storage root development under field and controlled greenhouse conditions. The transition of stem-derived roots towards starchstoring tuberous roots was analyzed via RNA sequencing (RNA-seq) and quantitative real-time PCR (qRT-PCR), as well as detailed histological studies. We identified pronounced auxin signaling events and the activation of secondary growth factors, as well as decreased GA signaling in the transition stages. This was accompanied by extensive regulation of cell wall biosynthesis genes, which probably is the result of the initial xylem expansion occurring during secondary growth within the root vasculature. After establishment of a vascular cambium, increased amounts of non-lignified parenchyma cells were observed. Starch storage was initiated in parenchyma cells and the corresponding samples displayed the transcriptional profiles of a starch storage crop. Interestingly, the activation of starch storage metabolism coincided with increased expression of the KNOX/BEL genes KNOTTED-LIKE FROM ARABIDOPSIS THALIANA (KNAT1), PENNYWISE (PNY) and POUND-FOOLISH (PNF), indicating their importance for proper xylem parenchyma function. Our study provides insight into the formation of the cassava vascular cambium and secondary growth, and provides a framework for subsequent genetic studies.

Planting material and growth conditions
Cassava stakes of genotype TME419 were planted in a field at IITA Ibadan, Nigeria towards the end of the rainy season. Root samples were taken from three individual stakes and frozen in liquid nitrogen at 30, 38, and 60 days after planting (dap). The samples were used for transcriptome analysis. Cassava stakes of genotype TME7 were grown in a greenhouse in Erlangen, Germany under a light regime of 12 h light and 12 h dark. Temperature was kept at a constant 30 °C and 60% relative humidity. Two nodal-derived root samples were taken from four stakes each at 22,26,30,34,38,42, and 60 dap. Approximately 5 mm root pieces of the primary bulking area at the proximal end of the root were stored in 70% ethanol for subsequent microscopy. Root tips were cut off and the root was frozen in liquid nitrogen. These samples were used for qRT-PCR.

Determination of soluble sugars, starch, and free amino acids
Soluble sugars, starch, and amino acids were measured as described previously (Obata et al., 2020).

Histology and microscopy
Histology and microscopy was performed as described previously (Mehdi et al., 2019).

RNA extraction, RNA sequencing, and qRT-PCR
Total RNA was extracted from TME419 roots by combining a modified cetyltrimethylammonium bromide (CTAB)-based extraction method (Li et al., 2008) with subsequent spin-column purification. Approximately 500 mg of sample material was ground in liquid nitrogen and mixed with 1 ml of pre-heated CTAB extraction buffer (2% CTAB, 2% PVP-40, 20 mM Tris-HCl, pH 8.0, 1.4 M NaCl, 20 mM EDTA). Samples were incubated at 65 °C for 15 min and centrifuged at 15 000 rpm at 4 °C for 5 min. The supernatant was transferred and mixed with an equal volume of cold chloroform:isoamyl alcohol (24:1) before centrifugation at 15 000 rpm for 10 min. The supernatant was mixed with 0.6 vol. of cold isopropanol and centrifuged at a maximum speed for 20 min. The pellet was washed with 70% ethanol, air-dried, and dissolved in nucleasefree water. After DNase I treatment, the resulting RNA was cleaned up using the kit RNA clean & concentrator™ (Zymo Research, USA) according to the manufacturer's instructions.
RNA samples were depleted of rRNA (Ribo-Zero rRNA Removal Kit Plant, Illumina) and sequenced with Illumina technology to obtain an average of 20 million paired-end reads. Raw files contained between 21 million and 60 million paired-end reads.
RNA extraction of TME7 roots was performed using the Spectrum Plant Total RNA Kit (Sigma-Aldrich, St. Louis, MO, USA). cDNA was generated from 0.5 µg of RNA using RevertAid H Minus Reverse Transcriptase as indicated by the manufacturer (Thermo Fisher Scientific, Waltham, MA, USA). The cDNA diluted was 1:10 and quantification of gene expression was examined using GoTaq® qPCR Master Mix (Promega, Madison, WI, USA). The assay was mixed in a 96-well plate and measured in an AriaMx Real-time PCR System (Agilent, Santa Clara, CA, USA). The results were analyzed using the 2 -ΔΔCt method (Livak and Schmittgen, 2001).

Results and discussion
Phenotypic and metabolic characterization of cassava root growth in the field Stakes of cassava genotype TME419 were planted in a field in Ibadan, Nigeria for up to 60 d to evaluate cassava storage root growth. The plants were carefully dug up and sampled at 30, 38, and 60 dap. Two types of roots were observed during the experiment, one of which progressed towards storage roots. These root types were subsequently called non-bulked roots (NBRs) and developing storage roots (DSRs). While the two root types were difficult to distinguish at 30 dap, they became distinct due to increased stiffness and darkening around 38 dap (Fig. 1A). This color change is likely to be the result of formation of the root periderm, a protective tissue layer that contains an enrichment of phenolic compounds (Campilho et al., 2020). The periderm is formed during secondary growth, making this root color change indicative for this process. Plants analyzed at 60 dap had clearly distintuishable tuberous roots that had already bulked considerably (Fig. 1A). Potential and visually distinguishable DSRs were analyzed for their soluble sugars, starch, proteins, and amino acids on a dry matter basis (Fig. 1B).
Soluble sugar concentrations increased linearly over time, with sucrose showing the highest abundance in comparison with glucose and fructose, indicating increased sucrose delivery and cleavage in storage roots. While there is a slight increase in starch concentration between 30 and 38 dap, an exponentially higher accumulation is observed between 38 and 60 dap. This is accompanied by a profound increase in free total amino acid concentrations. However, no elevation in protein concentration was observed. The high amino acid level at 60 dap was mainly due to accumulation of nitrogen-rich amino acids glutamine, glutamate, aspartate, and especially arginine, indicating increased delivery of nitrogen compounds from source tissues into storage roots (Fig. 1C). In contrast to many other root and tuber crops, cassava does not produce high abundant storage proteins (Vanderschuren et al., 2014), which could explain the high levels of free nitrogen-rich amino acids.
Together, the phenotypic and metabolic data suggested the onset of secondary root growth between 30 and 38 dap, and the onset of storage metabolism in the tuberous root between 38 and 60 dap.

Transcriptome analysis on cassava root types
To learn more about the underlying gene regulation of storage root formation, RNA-seq was performed on the different root samples. All replicates of each condition clustered together in a principal component anlysis (PCA) on log-transformed counts, demonstrating that no obvious sample outliers are present in the dataset ( Fig. 2A). The different time points and tissues separated along PC1 and PC2, respectively, with PC1 accounting for 64% and PC2 for 20% of variance. This indicated a higher influence of time than tissue on the cassava root transcriptome. The sampled roots separated into two different root types, DSRs or NBRs, for the respective time points (DSR30, DSR38, DSR60, NBR30, NBR38, and NBR60). The analysis indicated that the sampled root types were rather similar on their transcriptional level at 30 dap but became more different at later stages ( Fig. 2A). Only 90 of 3790 DEGs were specific for the time point 30 dap (Fig. 2B).
Most transcriptional changes were observed between the two root types at 38 dap, with >2000 genes being specifically regulated differently at this time point (Fig. 2B). This fitted the increasing visual distinction and onsetting secondary growth at 38 dap. Therefore, we focused our analysis on the remaining comparisons containing DSRs (DSR38/NBR38, DSR60/ NBR60, DSR38/DSR30, and DSR60/DSR38). The up-and down-regulated genes were used for KEGG pathway enrichment (see Table S2 at Dryad).
The pathway 'Sucrose and starch metabolism' (see Table  S3 at Dryad, sheet 1) was enriched in up-regulated genes of DSR60/NBR60 and DSR60/DSR38. Transcriptional up-regulation of starch biosynthesis during storage root development is expected considering the high starch content of cassava tuberous roots at 60 dap (Fig. 1B). The 'Amino sugar and nucleotide sugar metabolism' (Table S3 at Dryad, sheet 2) pathway was enriched in up-regulated genes of DSR38/ NBR38 and DSR60/NBR60. This pathway included the biosynthesis of diffferent ADP-, GDP-and UDP-sugars necessary for primary and secondary cell wall synthesis, which had their highest epression in storage roots at 38 dap. The enrichment at the time point 60 dap was mainly caused by some overlapping genes with the 'Sucrose and starch metabolism' pathway, as well as different α-glucosidases and chitinases. Fewer cell wall metabolism genes were observed in DSR60/NBR60.
Genes contained in 'Plant hormone signal transduction' (see Table S3 at Dryad, sheet 3) were solely enriched in DSR38/ NBR38, indicating that extensive developmental changes occurred around the 38 dap time point. Most notably, genes involved in auxin and cytokinin signaling were up-regulated in DSRs compared with NBRs. The significant enrichment of the 'Phagosome' (Table S3 at Dryad, sheet 4) in DSR38/ NBR38 was caused by different TUBULIN and RAC/RHO GTPase genes. Although these gene groups have many different functions, their regulation might reflect a high activity of the cytoskeleton during periods of cell division and cellular restructuring.

Regulation of starch and cell wall sugar biosynthesis
Starch and cell wall biosynthesis both compete for monosaccharide building blocks. This was also reflected in the differential expression of genes involved in either pathway Fig. 3. Overview of carbohydrate-related transcriptional changes between different root stages inferred from differentially expressed genes in the KEGG pathways 'Amino sugar and nucleotide sugar metabolism' and 'Sucrose and starch metabolism'. Heat maps show the log2FC of genes that are differentially expressed in at least one comparison. Columns from left to right are DSR38/NBR38, DSR60/NBR60, and DSR60/DSR38. Asterisks indicate a significant difference in expression, with |log2FC| ≥1 and FWER ≤0.05. Green rectangles highlight genes with differential expression. Light colored rectangles represent players with no differentially expressed genes in expression. Red rectangles indicate metabolites. (Fig. 3). Genes involved in cell wall formation were generally up-regulated in DSR38/NBR38, while genes involved in starch synthesis were generally up-regulated in DSR60/ NBR60. UDP-Glc is a key substrate in the biosynthesis of most secondary cell wall sugars, especially in the form of UDP-GlcA. This interconversion is catalyzed by the enzyme UDP-glucose dehydrogenase (UGDH) (Tenhaken and Thulke, 1996). One UGDH gene (Manes.14G084200) showed the most characteristic expression pattern in relation to the KEGG pathway analysis (Fig. 2C) by being strongly up-regulated in DSR38/ NBR38 and having next to no expression in DSR60. UGDH catalyzes the first non-reversible reaction in hemicellulose synthesis, and lack of the enzyme causes severe defects in secondary cell wall and pectin deposition as result of a low amount of UDP-GlcA (Reboul et al., 2011). UDP-GlcA serves as a substrate for UDP-Xyl synthesis through the enzymes UDP-xylose synthase (UXS) or UDP-apiose/UDP-xylose synthase (AXS) (Kuang et al., 2016). Two UXS genes (Manes.04G032900 and Manes.11G132500) were up-regulated at DSR38/NBR38 and either not up-regulated in DSR60/NBR60 or down-regulated in DSR60/DSR38. Xyl forms the backbone of xylans, the main polysaccharide in hemicellulose, whose biosynthesis is not yet fully understood (Pauly et al., 2013). However, cellulose synthase-like D (CSLD) enzymes take part in this process (Bernal et al., 2007). One CSLD gene (Manes.05G000900) was only up-regulated in DSR60/NBR60. These genes also play a role in correct meristem morphogenesis (Yang et al., 2016). Plants use UDP-GlcA to produce the main pectin subunit, UDP-GalA, via the enzyme UDP-glucuronic acid epimerase (GAE) (Gu and Bar-Peled, 2004). Pectin itself is in part synthesized from UDP-GalA through galacturonosyl transferase (GAUT) enzymes (Atmodjo et al., 2011). Five GAUT genes were differentially expressed and showed a similar trend to the other nucleotide-sugar biosynthesis genes (Manes.02G085400, Manes.03G122300, Manes.09G119300, Manes.14G088700, and Manes.15G076900). Pectins are part of the primary cell wall that exists in all cells. Generally speaking, nucleotide-sugar biosynthesis genes showed an up-regulation in DSR38/NBR38 and a down-regulation in DSR60/DSR38. This is especially true for genes directly using monosacharides UDP-Glc and fructose-6phosphate. This gene expression pattern was probably the result of different cellular compositions of the different samples, with 38 dap samples containing a higher ratio of cells displaying a secondary cell wall and samples of 60 dap containing more parenchyma cells.
One VACUOLAR INVERTASE (Manes.01G076500) and two CELL WALL INVERTASE (cwINV; Manes.03G049200 and Manes.11G025400) genes were down-regulated in a similar fashion. Reduction in cwINV expression fits the observation that cassava switches from an apoplasmic mode of transport in non-bulked roots to a symplasmic mode of transport in tuberous roots (Mehdi et al., 2019). Tuberous roots were shown to possess high sucrose synthase (SUS) activity (Mehdi et al., 2019). Indeed, the main SUS genes of tuberous roots (Manes.01G221900, Manes.03G044400, and Manes.16G090600) were highly expressed; however, their transcription was not significantly regulated between the different root samples. Three less expressed SUS genes, consisting of two SUS5 genes (Manes.01G123800 and Manes. 02G081500) and one SUS6 gene (Manes.14G107800), were down-regulated in DSR60/DSR38 and their proteins are only minor contibutors to the total SUS protein amount in tuberous roots (Mehdi et al., 2019). Similar to SUS, UDP-glucose pyrophosphorylase (UGPase), the enzyme that interconverts UDP-Glc and G-1-P, was highly expressed but not significantly differently regulated between the samples. Transcriptional up-regulation of SUS and UGPase was demonstrated in very early stages of potato tuber formation (Zrenner et al., 1993;Ferreira and Sonnewald, 2012), but could not be demonstrated for cassava tuberous roots in this study, indicating either an even earlier onset of SUS expression or post-transcriptional regulation.

Auxin signaling and vascular cambium establishment
Secondary growth in roots is accomplished through asymmetric cell division of stem cells. Previous studies on the morphology of cassava plants showed that adult storage roots continuously grow from a single vascular cambium ring, which resembles the secondary growth of model plants such as A. thaliana or poplar.
Other genes that were expressed along with PHB included important genes during vascular tissue development such as the receptor gene PHLOEM INTERCALATED WITH XYLEM (PXY) (Manes.14G037000) and its ligands CLAVATA3/ESR-RELATED 41 (CLE41) or CLE44 (Manes.12G052002 and Manes.13G045100). In A. thaliana roots, PXY is controlled through HD-ZIPIII activity and, together with CLE41/44, retains stem cell activity in the vascular cambium (Fisher and Turner, 2007;Hirakawa et al., 2008;Etchells and Turner, 2010;Smetana et al., 2019) through WUSCHEL-RELATED HOMEOBOX 4 (WOX4) and subsequent WOX14 expression. (Hirakawa et al., 2010;Etchells et al., 2013). However, no WOX4 gene was differentially expressed and no gene in the current annotation file was annotated as WOX14. Nevertheless, the closest AtWOX14 homolog (Manes.05G124300) was significantly up-regulated in DSR38/NBR38 (see Fig. S1 at Dryad). Further important genes in vascular cambium establishment, maintenance, and regulation of differentiation are KNOTTED-LIKE HOMEOBOX (KNOX) and BELLRINGER (BEL). Three KNAT1 genes in this dataset (Manes.05G184900, Manes.18G051700, and Manes.18G052901) were significantly up-regulated in DSR38/NBR38 and DSR60/NBR60, but showed much higher expression in DSR60 than in DSR38. KNOX proteins function as heterodimers together with BEL factors (Hay and Tsiantis, 2010). KNAT1 is known to interact with the functional paralogs PNY and PNF (Smith et al., 2004;Ragni et al., 2008), both of which have genes in cassava that showed the same expression pattern as KNAT1 (Manes.09G097500 and Manes.11G148000). These interaction pairs are known to repress the expression of boundary genes BLADE-ON-PETIOLE 1 (BOP1) or BOP2, as is the case here, with one BOP2 (Manes.08G056500) gene being down-regulated in DSRs (Woerlen et al., 2017). Despite this, the BOP2 gene was most highly down-regulated at DSR38/ NBR38 and the KNAT1 genes had their highest expression at 60 dap. Downstream factors of BOP2 showed a similar pattern to the boundary gene. These included a gene of the flowering gene APETALA1/AGAMOUS-like 7 (AP1/AGL7; Manes.01G103200). Another flowering time gene, AGL20 (Manes.05G041900), was down-regulated in a similar manner. Other crop species such as potatoes or onions utilize an adapted flowering pathway for storage organ initiation with genes of the FLOWERING LOCUS T (FT) serving as a mobile induction signal but also being expressed in the storage organ (Navarro et al., 2011;Lee et al., 2013;Abelenda et al., 2014). No expression of FT genes was found in any samples. The cassava storage root is a perennial organ that is not linked to reproduction. This implies a different regulation in comparison with vegetative reproduction organs-potatoes or onion bulbs. In addition, lack of flowering time genes has been linked to perennial growth (Melzer et al., 2008). KNAT1 was associated with the expression of the boundary gene LOB DOMAIN-CONTAINING PROTEIN (LBD) 4 . One gene was up-regulated in DSR38/NBR38 and DSR60/ NBR60 (Manes.06G173400). LBD4 takes part in keeping the procambium-phloem boundary together with PXY (Smit et al., 2020) and is able to promote the expression of PETAL LOSS (PTL) on the phloem side of the vasculature where it restricts expression of WOX4 to maintain organ boundaries . Three PTL genes were up-regulated in at least DSR38 versus NBR38 (Manes.10G055600), with two also being up-regulated at SR60/NBR60 and DSR60/ DSR38 (Manes.07G000300 and Manes.08G021500).
Overall, the results of the transcriptome analysis indicated that the DSRs underwent extensive cellular restructuring before the storage metabolism was activated. This can be derived from the strong transcriptional activation of genes encoding cell wall-modifying enzymes and hormone signaling components at stage 38 dap, as well as from the subsequent activation of storage metabolism between 38 and 60 dap. Although many of the developmental genes mentioned above were not significantly regulated within the comparisons of DSR samples, they showed differential expression between NBR and DSR samples, indicating that both tissues underwent a different program during their development. To gain a deeper insight into the underlying structural and transcriptional processes occurring in DSRs and to validate the field results independently, we performed a controlled greenhouse time-course experiment, covering the sampling time points of the RNA-seq experiment, as well as additional time points. The samples were used for histological analysis and qRT-PCR verification.

Morphological changes during root bulking
Stakes of cassava genotype TME7, a genotype closely related to TME419 (Bredeson et al., 2016), were planted in the greenhouse. Potential and developing storage roots emerging from the nodal areas of cassava stem pieces were sampled at seven time points (22, 26, 30, 34, 38, 42, and 60 dap) to obtain a broad line-up of DSRs at different ages and developmental stages. Cross-sections of the primary bulking site of every sampled root were taken, and the resulting sections were stained with toluidine blue solution in order to be able to group these roots according to their morphology, rather than their sampling time point (Fig. 5A-F).
Developing roots were sorted into five developmental stages based on their cellular composition and appearance. Stage one roots had a typical primary root morphology (Fig. 5A) and were observed in plants at 22 or 26 dap. They contained a starshaped central cylinder encased in an endodermis. Inside the central cylinder were alternating primary phloem and xylem cells with a procambium in between (indicated by a small layer of dividing cells).
The radial symmetry of the DSR was established between stages two and four (Fig. 5B-D). In stage two (Fig. 5B), the shape of the xylem and phloem tissue was less organized compared with stage one, indicating ongoing asymmetric cell division of the procambium that produces new phloem and xylem cells. The endodermis was pushed towards the outside, which led to a circular shape of the central cylinder. The cortex cells were still more or less intact. Stage three nodal roots showed a further increase in size of the central cylinder (Fig. 5C). The cambium was highly active and had a larger diameter than in the previous stage, albeit not yet being of radial shape. The newly formed cells were mostly lignified xylem cells, as indicated by the blue coloring of the newly formed cells towards the center. This fits the strong transcription activation of cell wall-related genes in field-grown plants at time point 38 dap. The phloem tissue did not increase much in comparison with younger roots. The endodermis was still intact, but the cortex was severely squashed towards the rhizodermis. The first vascular rays appeared but were hard to pinpoint at this stage. No starch was formed up to this point ( Fig. 5G-I). Stage four (Fig.  5D) DSRs were marked by dispersal of the endodermis and an annular shape of the inner xylem, as well as an overall radial shape of the root. The cortex was completely missing and the rhizodermis was displaced by a periderm at this stage. The root had not yet grown in diameter. However, the periderm gave it a visible brown coloring compared with the pale white of young roots. Furthermore, starch storage began in this phase of root development as indicated by Lugol staining (Fig. 5J). The stain was visible in the inner xylem cells adjacent to the vascular ray and the ray cells themselves, but most of the starch was still stored in phloem parenchyma cells. Similar to the field samples (Fig. 1A), a strong increase in starch is only measurable in established storage roots at ~60 dap (see Fig. S2A at Dryad). The color change between stage three and four was the main factor in root sampling for the transcriptome experiment, making this the likely developmental range of DSRs at 38 dap. Transcriptional changes at this time point were predominantly characterized by higher expression of secondary cell wall synthesis genes in tuberous roots. Secondary cell walls are solely deposited in differentiated sclerenchyma cells (Zhong et al., 2019). These include lignified xylem cells, which were strongly produced at this developmental stage (Fig. 5C, D). The number of meristematic cells increased, which explains the high expression of cambium genes, namely PXY (Etchells and Turner, 2010). Furthermore, the procambium developed into the ring-shaped vascular cambium during this time. This coincided with higher expression of auxin signal transduction genes and HD-ZIPIII. The more established storage roots (stage five; Fig. 5E, F) were easily identified on a macroscopic level due to an increase in lateral size. The vascular cambium mostly produced xylem parenchyma cells that were rapidly filled with starch, which is indicated by even the youngest cells being stained by the Lugol solution (Fig. 5K). This coincided with downregulation of secondary cell wall genes, up-regulation of starch biosynthesis genes, high amounts of starch (Figs 1B, 5K; see Fig. S2 at Dryad), as well as strong up-regulation of KNOX and BEL genes KNAT1, PNY, and PNF. In A. thaliana, KNAT1 is able to regulate the differentiation of xylem fibers in secondary meristems (Byrne et al., 2002;Mele et al., 2003;Liebsch et al., 2014) through regulation of GA signaling (Hay et al., 2002;Bolduc and Hake, 2009;Felipo-Benavent et al., 2018). In roots and hypocotyls of A. thaliana, KNAT1 promotes differentiation (Liebsch et al., 2014), while its overexpression in aerial parts of the plant leads to fewer lignified cells (Mele et al., 2003), a function that is similar in poplar trees (Groover et al., 2006). In other root and tuber crops, KNOX/BEL modules play a necessary role in the development of starchy storage organs. Sweet potatoes show similar expression levels of a BEL1 gene in late developmental stages compared with PNF in this study (Dong et al., 2019), which goes hand-in-hand with low GA signaling causing lower amounts of lignified cells to be produced (Singh et al., 2019). In potatoes, the KNOX gene POTATO HOMEOBOX 1 (POTH1) acts in a similar manner together with the BEL factor BEL5 (Rosin et al., 2003;Banerjee et al., 2006;Lin et al., 2013). Our study suggests that the low amount of lignified xylem cells was related to high KNAT1/PNF and/ or KNAT1/PNY activity in cassava tuberous root after establishment of the vascular cambium.

Independent qPCR validation of observed transcriptional changes
To validate and further elucidate the results of the RNA-seq analysis, a set of genes was chosen to represent the most important pathways and regulatory networks (Table 1; primer list  in Table S4 at Dryad) found to be regulated during storage root development in the transcriptome analysis. The expression was analyzed via qRT-PCR on the different DSR stages. The respective morphology and corresponding root stage were microscopically confirmed for every root.
The BEL factor genes PNY and PNF were strongly up-regulated in stage five DSRs, with a slight upward trend from stage one to four (Fig. 6). PNF expression was especially high, with a 40-fold up-regulation in the more established DSRs in comparison with stage one. A very similar trend is shown in the expression of KNAT1. The expression of these three genes correlated with the expression of starch genes PGM and PHS, but also with the production of a high number of xylem parenchyma cells (Fig. 5 E, F, K).
The boundary genes behaved as predicted. BOP2 expression showed a downward trend with especially low expression at stage five, which coincided with high expression of the KNAT1/BEL module. They also showed a significant negative correlation (Fig. 6, top left). LBD4, while not being up-regulated throughout root bulking, was expressed similarly to the transcriptome dataset at every stage.
PXY expression was detectable in every stage of root bulking. However, the expression dropped at stage three and went up again in later stages. This was highly correlated with the expression of PHB, which supports the hypothesis shown in Fig. 4  that PXY expression in cassava is dependent on an auxin signal inducing PHB expression, even though the change in the HD-ZIPIII gene is rather minute in comparison with PXY (Smetana et al., 2019;Zhang et al., 2019). This expression pattern was exactly contrary to the expression levels of GID1B (GIBBERELLIN INSENSITIVE DWARF 1B) and TPS as well as, to a lesser extent, GH3.1, PMI (PHOSPHOMANNOSE ISOMERASE), and BOP2. This coincided with deposition of lignified xylem cells by the vascular cambium between stage three and four (Fig. 5C, D) and with a comparably small up-regulation of the GA receptor gene GID1B, which might indicate a short GA signaling spike for the production of lignified cells. The starch biosynthesis genes PHS and PGM were strongly expressed at stage five in DSRs and have close to no expression at the other stages. This fits the histological analysis, where starch was mainly produced at stage five with negligible amounts in stage four and no starch in stages one to three ( Fig. 5G-L; see Fig. S2A at Dryad). Expression of these starch genes varied greatly within stage five roots but did correlate with the starch content of the particular roots (Fig. S2B, C at Dryad). AGL20 is strongly down-regulated between stage four and five, showing strong opposite behavior to the KNOX/ BEL and PGM/PHS genes. ANT, CLE41, and AUX1 increased in expression with advancing development. Similar to the transcriptomic results, SUS and UGPase did not exceed a |log2FC| of 1 between the DSR stages. The gene expression results obtained from independent samples with an independent method support the hypothesis formed through the transcriptome analysis, since the genes clustered in distinct groups (Fig. 6, top left) or otherwise resembled the transcriptional profile shown in Figs 3 and 4. GID1B, TPS, BOP2, and GH3.1 formed a cluster as expected and were more highly expressed in early stages of nodal roots (Fig. 6, cluster 1). The proposed auxin signaling genes (AUX1, PHB, and PXY) were tightly correlated with each other, as well as with ANT (Fig. 6, cluster 2). Starch biosynthesis genes were only strongly expressed in stage five roots together with the KNOX/BEL genes (Fig. 6, cluster 6), in accordance with the transcriptome dataset.

Conclusion
This study investigated the transcriptional and morphological changes occurring during early storage root development. We show that the lateral size increase in storage roots was accomplished through establishment of a singular vascular cambium and subsequent deposition of non-lignified parenchyma cells. Early storage root development was marked by production of lignified cells in the central cylinder. This coincided with the strong expression of important secondary cell wall biosynthesis genes, namely UGDH, GAE, UXS, and GAUT. The same genes were down-regulated in later stages where starch biosynthesis is the dominantly active pathway. Starch biosynthesis began only after the establishment of the vascular cambium, after which mainly xylem parenchyma cells were produced. We demonstrate that this was accompanied by strong transcriptional up-regulation of plastid-localized starch biosynthesis genes such as PGM, AGPase, GBSS, SS, and BE. The vascular cambium itself was formed through a transcriptionally strong auxin response correlated with expression of a HD-ZIPIII gene and subsequent expression of PXY, CLE41/44, and WOX14. The formation of parenchyma rather than sclerenchyma cells was accompanied by high expression of the KNOX/BEL genes KNAT1, PNY, and PNF, highlighting their importance for xylem parenchyma function.
While this study provides a groundwork for a lesser known crop species, it also raises many questions to be answered. For instance, how will direct genetic alteration of different developmental regulators alter storage root expansion and storage parenchyma formation? Understanding the role of KNOX/ BEL factors and the role of different hormones in the suppression of secondary cell wall biosynthesis in storage parenchyma cells should yield useful information for further advances in cassava.