Cassava storage roots are among the most important root crops worldwide and represent one of the most consumed staple foods in Sub-Saharan Africa. The vegetatively propagated tropical shrub can form many starchy tuberous roots from its stem. These storage roots are formed through the activation of secondary root growth processes. However, the underlying genetic regulation of storage root development is largely unknown. Here we report on distinct structural and transcriptional changes occurring during the early phases of storage root development. A pronounced increase in auxin-related transcripts and the transcriptional activation of secondary growth factors, as well as a decrease in gibberellin-related transcripts was observed during the early stages of secondary root growth. This was accompanied by increased cell wall biosynthesis, increased most notably during the initial xylem expansion within the root vasculature. Starch storage metabolism was activated only after the formation of the vascular cambium. The formation of non-lignified xylem parenchyma cells and the activation of starch storage metabolism coincided with increased expression of the KNOX/BEL genes KNAT1, PENNYWISE and POUND-FOOLISH, indicating their importance for proper xylem parenchyma function.

Planting material and growth conditions

Cassava stem sticks 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 sticks and frozen in liquid nitrogen at 30 dap, 38 dap and 60 dap. The samples were used for transcriptome analysis. Cassava stem sticks of genotype TME7 were grown in a green house in Erlangen, Germany under a light regime of 12 h light and 12 h dark. Temperature was kept at a constant of 30°C and 60% relative humidity. Two nodal- and two cambium-derived root samples from the basal end of the stick were taken from four sticks 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% EtOH 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 CTAB-based extraction method (Li et al., 2008) with subsequent spin-column purification. Approximately 500mg of sample material was grinded in liquid nitrogen and mixed with pre-heated 1 mL 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 15000 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 15000 rpm for 10 min. The supernatant was mixed with 0.6 volume of cold isopropanol and centrifuged at maximum speed for 20 min. 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, USA) according to manufacturer's instructions.

RNA samples were depleted of ribosomal RNA (Ribo-Zero rRNA Removal Kit Plant, Illumina) and sequenced with Illumina technology to obtain an average of 20 million paired-end reads. Raw files containing 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. Luis, MO, USA). cDNA was generated from 0.5 µg RNA using the RevertAid H Minus Reverse Transcriptase as indicated by the manufacturer (Thermo Fisher Scientific, Waltham, MA, USA). The cDNA was 1:10 diluted and quantification of gene expression was examined using GoTaq® qPCR Master Mix (Promega, Madison, USA). The assay was mixed in a 96-well plate and measured in an AriaMx Real-time PCR System (Agilent, Santa Clara, USA). The results were analyzed using the 2^-ΔΔCt method (Livak and Schmittgen, 2001).

Read trimming and mapping

FastQ files containing the raw sequencing reads were quality checked using FastQC (v. 0.11.5; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and MultiQC (v. 1.8; https://multiqc.info/). Adapter and quality trimming was performed in two steps utilizing the k-mer trimming tool BBduk (v. 38.96; https://sourceforge.net/projects/bbmap/) with its provided adapter sequences. A k-mer length of 21 was set allowing a minimum k-mer length of 11 and two mismatches. Reads < 35 nucleotides or an average quality < 20 were excised, as well as individual bases below a quality of 20 at the ends of the read. The resulting FastQ files were mapped to the M. esculenta genome (v.7.1; https://genome.jgi.doe.gov/portal/pages/dynamicOrganismDownload.jsf?organism=Mesculenta) in two passes using STAR (v.2.5.0a; Dobin et al. (2012); https://github.com/alexdobin/STAR).  The resulting BAM files were indexed and deduplicated employing samtools (v.1.7; Li et al. (2009) ; http://www.htslib.org/). Read counting was performed using the program featureCounts (v.1.5.0; Liao et al. (2013); http://bioinf.wehi.edu.au/featureCounts/). Only primary reads were counted. Trimmed, mapped and deduplicated read counts are available in table S1. All aforementioned programs were used under Linux (Ubuntu v. 18.04 LTS).

Data analyses

Log2 fold-change (log2FC) and its standard error were estimated in R (v. 3.6.2) utilizing the Bioconductor package DESeq2 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html; Love et al. (2014)) on individual pairs. Wald’s test was used to calculate p-values between pairs, which were adjusted after Bonferroni’s family wise error rate (FWER). Genes with |log2FC| ≥ 1 and FWER ≤ 0.05 were accepted as differentially expressed genes (DEGs). Enrichment analysis were conducted with a one-sided Fisher’s exact test using the Bioconductor package clusterProfiler. (https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html; Yu et al. (2012)). Enrichments with FWER ≤ 0.05 were accepted as significant. Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology (KO) terms, cassava and tale cress identifiers were taken from an annotation file published with the genome. Pathway and regulatory networks were constructed through publication- and database mining (STRING [https://string-db.org/], BioGRID [https://thebiogrid.org/] and TAIR [https://www.arabidopsis.org/]). In the text, cassava genes were described by their best A. thaliana hit based on BLASTP similarity.