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

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.

Methods
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.

木薯贮藏根是全球最重要的块根作物之一，亦是撒哈拉以南非洲地区最主要的主食作物之一。这种以营养繁殖的热带灌木可从茎部生成大量淀粉质块根。此类贮藏根的形成依赖于次生根系生长过程的激活，但目前学界对其发育背后的遗传调控机制仍知之甚少。本研究报道了贮藏根发育早期阶段发生的独特结构与转录组变化。在次生根系生长的早期阶段，我们观察到生长素（auxin）相关转录本显著上调，次生生长相关因子的转录激活，以及赤霉素（gibberellin）相关转录本的下调；同时伴随细胞壁生物合成增强，这一增强效应在根维管组织的初始木质部（xylem）扩张阶段尤为显著。淀粉贮藏代谢仅在维管形成层（vascular cambium）形成后才被激活。非木质化木质部薄壁细胞的形成与淀粉贮藏代谢的激活，与KNOX/BEL基因（KNOX/BEL genes）的KNAT1、PENNYWISE及POUND-FOOLISH的表达上调同时发生，表明这些基因对木质部薄壁细胞的正常功能具有重要作用。

方法

## 实验材料与生长条件
将基因型为TME419的木薯茎秆插条种植于尼日利亚国际热带农业研究所（IITA）伊巴丹校区的田间，种植时间为雨季末期。分别于移栽后天数（days after planting, dap）30、38和60时，从3个独立茎秆插条采集根样本，经液氮速冻后用于转录组分析。
将基因型为TME7的木薯茎秆插条种植于德国埃尔朗根的温室中，光照周期设置为12小时光照/12小时黑暗，温度恒定为30℃，相对湿度维持在60%。分别于22、26、30、34、38、42和60 dap时，从4个茎秆插条的基部采集2个节部来源根样本与2个形成层来源根样本。采集根近端区域约5mm的膨大初生区域样本，保存于70%乙醇中用于后续显微观察；切取根尖后将剩余根组织经液氮速冻，此类样本用于定量实时荧光定量PCR（quantitative real-time PCR, qRT-PCR）分析。

## 可溶性糖、淀粉与游离氨基酸含量测定
可溶性糖、淀粉及氨基酸的含量测定参照此前报道的方法（Obata等，2020）。

## 组织学与显微观察
组织学与显微观察操作参照此前报道的方法（Mehdi等，2019）。

## RNA提取、RNA测序与qRT-PCR分析
针对TME419根样本的总RNA提取，采用改良十六烷基三甲基溴化铵（hexadecyltrimethylammonium bromide, CTAB）提取法（Li等，2008）结合后续离心柱纯化的组合方案：称取约500mg样本于液氮中研磨粉碎，加入1mL预热的CTAB提取缓冲液（2% CTAB、2% PVP-40、20 mM Tris–HCl pH 8.0、1.4 M NaCl、20 mM EDTA），65℃孵育15分钟后，4℃下15000 rpm离心5分钟。将上清液转移并加入等体积预冷的氯仿:异戊醇（24:1）混合，随后15000 rpm离心10分钟。取上清液加入0.6倍体积的预冷异丙醇，最大转速离心20分钟。沉淀用70%乙醇洗涤，风干后溶于无核酸酶水中。经DNase I处理后，使用RNA clean & concentrator™试剂盒（"Zymo Research", 美国）按照制造商说明书进行RNA纯化。
RNA样本经Ribo-Zero rRNA Removal Kit Plant（Illumina）去除核糖体RNA后，采用Illumina测序平台进行测序，平均每个样本获得2000万条双端读段。原始测序文件包含2100万至6000万条双端读段。
针对TME7根样本的RNA提取采用Spectrum植物总RNA提取试剂盒（"Sigma-Aldrich", 美国密苏里州圣路易斯市）。取0.5 μg RNA作为模板，使用RevertAid H Minus逆转录酶（"Thermo Fisher Scientific", 美国马萨诸塞州沃尔瑟姆市）按照制造商说明书合成cDNA。将cDNA按1:10比例稀释后，采用GoTaq® qPCR Master Mix（"Promega", 美国麦迪逊市）进行基因表达定量分析。反应体系置于96孔板中，使用AriaMx实时荧光定量PCR系统（"Agilent", 美国圣克拉拉市）进行检测。结果采用2-ΔΔCt法进行分析（Livak与Schmittgen，2001）。

## 读段修剪与比对
原始测序读段的FastQ文件质量使用FastQC（v.0.11.5; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/）与MultiQC（v.1.8; https://multiqc.info/）进行评估。接头与质量修剪分为两步进行，使用k-mer修剪工具BBduk（v.38.96; https://sourceforge.net/projects/bbmap/）并采用其自带的接头序列。设置k-mer长度为21，允许最小k-mer长度为11，允许2个错配。移除长度小于35核苷酸或平均质量值低于20的读段，同时切除读段末端质量值低于20的单个碱基。将处理后的FastQ文件分两轮比对至木薯（Manihot esculenta, M. esculenta）参考基因组（v.7.1; https://genome.jgi.doe.gov/portal/pages/dynamicOrganismDownload.jsf?organism=Mesculenta），比对工具采用STAR（v.2.5.0a; Dobin等，2012; https://github.com/alexdobin/STAR）。生成的BAM文件使用samtools（v.1.7; Li等，2009; http://www.htslib.org/）进行索引与去重复。读段计数采用featureCounts程序（v.1.5.0; Liao等，2013; http://bioinf.wehi.edu.au/featureCounts/）完成，仅统计初级读段。经过修剪、比对与去重复的读段计数结果可见表S1。所有上述软件均运行于Linux系统（Ubuntu v.18.04 LTS）。

## 数据分析
在R（v.3.6.2）环境中，利用Bioconductor包DESeq2（https://bioconductor.org/packages/release/bioc/html/DESeq2.html; Love等，2014）对各组样本对估计log2倍变化（log2 fold-change, log2FC）及其标准误。采用Wald检验计算样本对之间的p值，并经Bonferroni法校正家族式错误率（family wise error rate, FWER）。将|log2FC|≥1且FWER≤0.05的基因认定为差异表达基因（differentially expressed genes, DEGs）。富集分析采用单侧Fisher精确检验，通过Bioconductor包clusterProfiler（https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html; Yu等，2012）完成，将FWER≤0.05的富集结果认定为显著富集。京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes, KEGG）直系同源（orthology, KO）术语、木薯与拟南芥（Arabidopsis thaliana, A. thaliana）基因标识符均取自该基因组发布时配套的注释文件。通路与调控网络通过文献与数据库挖掘构建（STRING [https://string-db.org/], BioGRID [https://thebiogrid.org/]和TAIR [https://www.arabidopsis.org/]）。正文中木薯基因以其基于BLASTP相似性比对得到的最优拟南芥同源基因进行描述。