https://orcid.org/0000-0002-6953-7995
Abstract
Cassava (Manihot esculenta) is an important starchy root crop that provides food for millions of people worldwide, but little is known about the regulation of the development of its tuberous root at the multi-omics level. In this study, the transcriptome, proteome, and metabolome were examined in parallel at seven time-points during the development of the tuberous root from the early to late stages of its growth. Overall, highly dynamic and stage-specific changes in the expression of genes/proteins were observed during development. Cell wall and auxin genes, which were regulated exclusively at the transcriptomic level, mainly functioned during the early stages. Starch biosynthesis, which was controlled at both the transcriptomic and proteomic levels, was mainly activated in the early stages and was greatly restricted during the late stages. Two main branches of lignin biosynthesis, coniferyl alcohol and sinapyl alcohol, also functioned during the early stages of development at both the transcriptomic and proteomic levels. Metabolomic analysis further supported the stage-specific roles of particular genes/proteins. Metabolites related to lignin and flavonoid biosynthesis showed high abundance during the early stages, those related to lipids exhibited high abundance at both the early and middle stages, while those related to amino acids were highly accumulated during the late stages. Our findings provide a comprehensive resource for broadening our understanding of tuberous root development and will facilitate future genetic improvement of cassava.
Introduction
Cassava (Manihot esculenta) is an important tropical/sub-tropical starchy root crop, providing staple food for >750 million people worldwide (Yang et al., 2011). The main economic value of cassava is in its tuberous roots, which can be made up of >80% starch on a dry weight basis (Li et al., 2010), and hence much research has focused on its root development. Cassava is usually propagated by stem cuttings, from the base of which fibrous roots are regenerated and then initiate into root tubers from 6–16 weeks after planting, with the main period being at 3–4 months (Lowe et al., 1982). The root tubers, which became visible to the naked eye around the 8th week, increase in size via secondary growth with an increase in width but not in length, and this is accompanied by a decrease in sugar content and an increase in starch content (Lowe et al., 1982). Following this root bulking, carbohydrates synthesized in the leaves are translocated into the roots mainly during the period 80–300 days after planting, although this can vary among different genotypes and/or environmental conditions (El-Sharkawy, 2004). Corresponding to the pattern of development, the dry matter of the roots increases slowly during the first 2 months after planting, then increases rapidly from 3–10 months, after which it is maintained at a high level for further 2 months (El-Sharkawy, 2004). The phenotypic changes that occur during root development might result from stage-dependent regulation of metabolic pathways, as has been reported in other species (Casas-Vila et al., 2017; Slade et al., 2017).
Many metabolic pathways have been shown to be associated with root development in cassava. The most significant changes are in carbohydrate and energy metabolism, followed by secondary metabolism, and amino acid and lipid metabolism (Li et al., 2010). Yang et al. (2011) reported similar results, but suggested that glycolysis/gluconeogenesis is the most important pathway because the component exchanges with other pathways are strongly dependent on its existence. Cell wall biosynthesis, which plays important roles in root elongation and expansion, also shows significant changes during the development of tuberous roots (Yang et al., 2011; Wang et al., 2016b), and corresponding alterations are found in the metabolism of lignin, which functions in the reinforcement of cell walls (Gray et al., 2012). In addition, transcription factors and hormone-related genes also play essential roles in the initiation, formation, and development of tuberous roots (Yang et al., 2011; Sojikul et al., 2015; Dong et al., 2019). Thus, the development of tuberous roots in cassava is a very complex process that involves changes in many metabolic pathways, and there is a need for the underlying regulatory mechanisms to be systematically studied.
Transcriptional and post-transcriptional regulation both play crucial roles in root development and starch metabolism (Zhang et al., 2009; D’haeseleer et al., 2011; Dong et al., 2019). MtNAC1 is a key transcription factor in Medicago truncatula that is involved in lateral root formation via both transcriptional and post-transcriptional regulation (D’haeseleer et al., 2011). In sweet potato (Ipomoea batatas), many genes related to the initiation and development of storage roots are also regulated at the transcriptional and post-transcriptional levels (Dong et al., 2019). In starch synthesis in wheat (Triticum aestivum), SS, SBE, and DBE are regulated transcriptionally, whereas GBSS-I is regulated post-transcriptionally (Wang et al., 2014). In cassava, post-transcriptional silencing of SBE2 results in very high contents of amylose in roots (Zhou et al., 2019).
A number of studies have been performed using high-throughput technologies to investigate the regulation of cassava root development. At the transcriptomic level, Yang et al. (2011) identified thousands of differentially expressed genes related to the initiation and formation of cassava storage roots, Li et al. (2010) developed a cDNA microarray to analyse gene expression profiling across four key growth stages, whilst Sojikul et al. (2015) determined that hormones and their interactions play essential roles in storage root initiation. At the proteomic level, Sheffield et al. (2006) identified 237 proteins differentially accumulated in fibrous and tuberous roots, Mitprasat et al. (2011) suggested that possible metabolic switches in leaves might regulate/trigger storage root initiation and growth, and Wang et al. (2016b) detected 154 differentially expressed proteins during starch accumulation and root tuberization. However, these studies were each performed either at the transcriptomic or proteomic level and could not differentiate transcriptional and post-transcriptional regulation in root development. At the metabolomic level, no studies have been performed that examine tuberous root development in cassava, although there have been reports on post-harvest physiological deterioration (Uarrota et al., 2014), descriptions of biochemical diversity (Drapal et al., 2019; Obata et al., 2020; Price et al., 2020), and quantification protocols (Rosado-Souza et al., 2019).
Multi-omics integrative analysis is a promising tool to systematically examine complex physiological processes (Garg et al., 2019; Hu et al., 2019; Li et al., 2019). In this study, we investigated the global abundances of mRNAs, proteins, and metabolites in parallel at seven developmental time-points during tuberous root development in cassava, and found highly dynamic and stage-specific changes. We then used a multi-omics integrative analysis to explore the stage-specific roles and the regulatory pathways of the genes/proteins. Our results provide a comprehensive resource to expand our understanding of tuberous root development and will facilitate future genetic improvement of cassava.
Materials and methods
Plant growth and sample collection
As described previously (Ding et al., 2019), stems of cassava (Manihot esculenta cv. ‘SC205’) were cut into small sections of ~15 cm length with 2–3 buds, and grown under field conditions at the experimental farm of the Chinese Academy of Tropical Agricultural Sciences, DanZhou (109.50°E, 19.51°N), with normal crop management practices. The daily mean temperature ranged from 24–35 °C. Seven experimental blocks were established for sampling at each of seven time-points. Each block consisted of four rows with each row having seven plants. Only the five plants growing in the middle of the two central rows were used for sampling, which took place in the morning (~09.00 h). The seven time-points (S1–S7) were 100, 140, 180, 220, 260, 300, and 340 d after planting (DAP). These were considered as the early (S1–S3), middle (S4), and late (S5–S7) stages of growth (Li et al., 2010). When the plants were small (S1, S2) six replicate plants were usually sampled, and five were sampled when the plants were bigger (S3–S7). The plants were dug up and one typical main root per plant was selected for analysis.. The root was sectioned from the middle, cut into slices ~3 mm thick, and 5–6 of the slices were frozen immediately in liquid nitrogen and kept at −80 °C until use.
Root diameter was measured for the 5–6 replicate plants. Starch and soluble sugar contents were measured as described by Wang et al. (2016b) with three biological replicates. The characteristics of starch granules in the roots were examined using a Hitachi S-3000N SEM with at least two biological replicates.
Transcriptomic analysis
RNA sequencing (RNA-seq) and library construction were performed by the Annoroad Gene Technology Corporation (Beijing, China). Briefly, ~3 μg RNA from each sample was extracted using a TruSeq RNA Sample Prep Kit (Illumina), according to the manufacturer’s instructions. RNA-seq was conducted on an Illumina Hiseq 4000 to produce reads 2 × 150 bp in length. Samples from three biological replicates were used.
Sequence quality was examined using the FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), as described previously (Ding et al., 2016). Low-quality reads and adaptor contaminations were removed using FASTX-toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html). Clean reads were aligned to the cassava reference genome (version 6.1) obtained from the Phytozome database (https://phytozome.jgi.doe.gov/pz/portal.html) using TopHat v2.0.10 (Kim et al., 2013), and then assembled with a reference genome-based strategy using Cuflinks v2.1.1 (Trapnell et al., 2012). Differentially expressed genes (DEGs) were determined using DESeq2 (Love et al., 2014) with settings of log2|fold-change| >1 and false discovery rate <0.05. Gene expression was normalized as fragments per kilobase per million mapped reads (FPKM).
Proteomic analysis
The proteome was determined by Shanghai Luming Biological Technology Co., Ltd, using iTRAQ as previously described (Ding et al., 2019). Briefly, total protein was extracted from each sample, and the concentration was determined using a Bio-Rad Protein Assay kit with bovine serum albumin (BSA) as the standard. After trypsin digestion, the peptides were dried using vacuum centrifugation and labeled with 8-plex iTRAQ reagents (Applied Biosystems) following the manufacturer’s protocol.
For LC-MS/MS assays, the iTRAQ-labeled peptide mix was loaded onto a C18 column (100 µm × 30 mm, 3 µm, 150 Å) and then washed with Nano-RPLC buffer A (2% ACN and 0.1% FA) for 10 min at 2 μl min–1. Acquisition of the MS data was performed on using a Triple TOF 5600 system (Sciex) with parameters described previously (Ding et al., 2019).
The MS/MS data generated were searched against the cassava reference database for protein identification using ProteinPilot v5.0 (Sciex). The Proteomics System Performance Evaluation Pipeline (PSPEP) was applied to calculate the false discovery rate based on an automatic decoy database search. Unique peptides with confidence higher than 95% were kept for iTRAQ quantification, and proteins with an unused value >1.3 were retained for further analysis. Differences in protein abundance were compared using Student’s t-test with Benjamini–Hochberg correction. Differentially expressed proteins (DEPs) were determined using a fold-change >1.2 (or <0.83) and q-value <0.05. Samples from two biological replicates were used.
Metabolic analysis
The untargeted metabolome was determined by Wuhan Metware Biotechnology Co., Ltd, as described previously (Chen et al., 2013). Briefly, 100 mg of powdered sample was extracted overnight with 1.2 ml 70% aqueous methanol at 4 °C. After centrifugation for 10 min at 10 000 g, the extracts were absorbed and filtrated. A UPLC-ESI-MS/MS system (UPLC, Shim-pack UFLC SHIMADZU CBM30A; MS, Applied Biosystems 6500 Q TRAP) was used to analyse the extracts with following conditions: UPLC column, Waters ACQUITY UPLC HSS T3 C18 (1.8 µm, 2.1 mm × 10 cm). The mobile phase consisted of solvent A (pure water with 0.04% acetic acid) and B (acetonitrile with 0.04% acetic acid). Sample measurement was conducted with a gradient program employing start conditions of 95% A and 5% B. A linear gradient to 5% A and 95% B within 10 min was programmed, and a composition of 5% A and 95% B was kept for 1 min. Subsequently, the composition was adjusted to 95% A and 5% B within 0.1 min and kept for 2.9 min. The injection volume was set to 2 μl and the column oven was 40 °C. The effluents were connected to an electrospray ionization (ESI)-triple quadrupole-linear ion trap (Q TRAP)-MS.
Acquisition of the linear ion trap (LIT) and triple quadrupole (QQQ) scans was performed on a triple quadrupole-linear ion trap mass spectrometer API 6500 Q TRAP LC/MS/MS system, equipped with a Turbo Ion-Spray interface (operating in positive and negative ion mode) and manipulated using Analyst 1.6.3 (Sciex). The parameters of the ESI source operation were set as follows: ion source, turbo spray; source temperature, 550 °C; ion spray voltage, 5500 V for positive ion mode and –4500 V for negative ion mode; ion source gas I, gas II, and curtain gas were set to 50, 60, and 30 psi, respectively; the collision gas was high. Instrument tuning and mass calibration were conducted in QQQ and LIT modes using 10 μmol l–1 and 100 μmol l–1 polypropylene glycol solutions, respectively. The QQQ scan was acquired during the MRM experiment setting the collision gas (nitrogen) equal to 5 psi.
Metabolites were identified using information in public metabolite databases and the Metware database of the Metware Biotechnology Co. All identified metabolites were subjected to principle component analysis (PCA) and orthogonal partial least squares discriminate analysis (OPLS-DA), and significant differences were determined by setting the variable importance in projection (VIP) as ≥1 and log2 |fold-change| ≥1. Samples from three biological replicates were used.
Multi-omics integrative analysis
The abundances of mRNAs, proteins, and metabolites were min–max normalized between –1 and 1 among the different time-points (Becker et al., 2018). The abundance patterns were determined using the standard procedure of the WGCNA package in R (Langfelder and Horvath, 2008; Ding et al., 2015) and subsequently visualized using the R package ‘pheatmap’. mRNAs, proteins, and metabolites that were non-assigned into any groups were excluded for further analysis.
To interpret the biological functions of the DEGs and DEPs, cassava genes were functionally annotated and assigned into different hierarchical categories using the MapMan classification system (Thimm et al., 2004). The significances of the enriched categories were using by Fisher’s exact test, as previously described (Ding et al., 2016; Fu et al., 2016). Pearson’s correlation test was applied to calculate the correlation coeficients between the DEPs and their corresponding mRNAs. WGCNA was also used to determine the associations between mRNAs/proteins and metabolites, and Cytoscape (Su et al., 2014) was used for network visualization.
Results
Phenotypic changes during root development
To investigate phenotypic changes during tuberous root development in cassava, roots were harvested every 40 d from 100 DAP to 340 DAP for a total of seven time-points, referred to as S1–S7 and corresponding to the early stage (S1–S3), middle stage (S4), and late stage (S5–S7) of growth , (Fig. 1A). The root diameter steadily increased from S1 to S4, followed by a sharp increase at S5, and then a further slight increase to S7 (Fig. 1B). The starch content initially decreased from S1 to S2 but then steadily increased from S4 onwards, while the soluble sugar content decreased significantly from S1 to S5 before increasing at S6 (Fig. 1C). SEM showed that starch granules were clearly visible at S2 (Fig. 1D) and S3 (Fig. 1E), and the numbers dramatically increased at S4 (Fig. 1F) and S6 (Fig. 1G), indicating the importance of stage S4 for starch accumulation.
Phenotypic variation during tuberous root development in cassava. (A) Sections of roots sampled at seven time-points from 100–340 d after planting: S1, 100; S2, 140: S3, 180: S4, 220; S5, 260; S6, 300; and S7, 340. S1–S3 represents the early stage of root development, S4 represents the middle stage, and S5–S7 represents the late stage. The scale bar is 2 cm. (B) Root diameters and (C) contents of starch and soluble sugars at the different developmental stages. Data are means (±SD) from at least three biological replicates. Different letters indicate significant differences between means as determined using ANOVA followed by Duncan’s multiple range test (P<0.05). (D–G) SEM images of the middle parts of tuberous roots collected at (D) S2, (E) S3, (F) S4, and (G) S6. Scale bars are 10 μm.
Root samples collected from the seven time-points were subjected to transcriptomic, proteomic, and metabolomic analyses in parallel. Overall, the sample clustering revealed highly dynamic and coordinated changes in the abundances of mRNAs, proteins, and metabolites during root development (Fig. 2A, Supplementary Fig. S1 at JXB online), for example S1–S3 were closely clustered and S5–S7 were closely clustered, while S4 showed an intermediate relationship.
Dynamic, coordinated, and stage-specific changes in mRNAs, proteins, and metabolites during tuberous root development in cassava. (A) Clustering of samples for significantly changed mRNAs, proteins, and metabolites at the seven time-points during development (see Fig. 1). (B) Profiles of genes identified as being differentially expressed during development exclusively at the transcriptomic level. A total of 12 groups were determined according to the pattern of expression (gM1–gM12). (C) Profiles of genes identified as being differentially expressed at during development at both the transcriptomic and proteomic levels. Four groups were determined (gpM1–gpM4). (D) Profiles of genes identified as being differentially expressed during development exclusively at the proteomic level. Two groups were determined (pM1–pM2). (E) Changes in abundance of metabolites during development. Six groups were determined (mM1–mM6). Samples/groups are identified as follows: g, transcriptomic (mRNA) level; p, proteomic level; and m. metabolomic level.
Transcriptomic and proteomic profiling during root development
In total, ~853.5 million clean reads were generated by RNA-seq, of which 83.8% were aligned to the cassava reference genome and used to generate an expression matrix across samples (Supplementary Tables S1, S2). A threshold of FPKM>1 was arbitrarily selected to identify expressed genes among the samples. A total of 22 463 genes were expressed across all seven time-points, which was ~68% of the annotated genes in the reference genome. Interestingly, the number of expressed genes gradually decreased as the roots developed from S1 (20 250) to S7 (18 166).
A total of 303 587 mass spectra were obtained from the root samples using an iTRAQ approach. By removing low-scoring spectra and searching against cassava reference proteins, a total of 13 209 unique peptides were generated and these were further grouped into 3970 proteins. After setting the unused value at >1.3 and peptides at ≥1, a total of 2657 non-redundant proteins were confidently determined and quantified (Supplementary Table S3).
Integrative transcriptomic and proteomic analysis
An integrative analysis was performed between the transcriptomic and proteomic data, and 8556 DEGs and 564 DEPs were exclusively identified, with 805 being commonly identified between the two sets of data. We clustered these DEGs and DEPs according to their expression patterns, and conducted functional enrichment analysis.
DEGs exclusively identified at the transcriptomic level
A total of 12 groups of genes (gM1–gM12) were identified that were significantly differentially expressed among the seven stages at the mRNA level but not at the protein level (Fig. 2B, Supplementary Table S4).
Genes from groups gM1–4 were mainly expressed at the early stages of tuberous root development. For example, gM1 had highest expression at S1 and then steadily decreased until S7. The genes of this group were significantly enriched in auxin and pathways related to the cell wall (Fig. 3A, Supplementary Table S5). Many genes included in this group were related to cellulose synthesis (CESA1, CESA3, CESA6, and CSLC5), cell wall proteins (FLA1, FLA4, RLP29, and RGP3), synthesis of cell wall precursors (GAE3, GAE6, RHM1, and UXS2), and cell wall modification (EXPA1, EXPA8, XTH5, and XTH8). Several genes related to auxin biosynthesis (ILL3 and ILL6), signal transduction (PIN1, ABP1, and TIR1), and auxin responses (DFL1, AIR9, RRT1, SAUR12, and SAUR14) were also included (Supplementary Fig. S2). These results strongly indicated a coordinated contribution of both cell wall and auxin genes to early development in the tuberous roots. Genes from gM2–4 generally had highest expression at S3, although gM2 showed relatively high expression at S1 and S2 while gM4 exhibited relatively high expression at S4. Genes in gM3 and gM4 were enriched in abscisic acid (ABA) and abiotic stress, and in ethylene metabolism, respectively, whereas no enrichment was detected for gM2 (Fig. 3A). Several genes involved in ABA biosynthesis (NCED3 and NCED9), ABA signaling (ABF3 and ABI1), ABA responses (HAI3 and HVA22E), and heat-shock proteins (HSP20, HSP17.4, and HSP90.1) were included in gM3 (Supplementary Table S4). A few genes participating in ethylene biosynthesis (ACS6) and ethylene signaling transduction (ERF5-ERF7 and ERF105) were found in gM4 (Supplementary Table S4). These enrichments suggested that hormone signaling has important roles in early the development of the tuberous roots.
Functional enrichment of genes/proteins and their association with metabolites during development of tuberous roots in cassava. (A) Enrichment of functional categories in groups of genes/proteins identified at the transcriptomic and/or proteomic levels. (B) Associations between the groups of genes/proteins (rows) and the groups of metabolites (columns). Each cell contains the corresponding Pearson’s correlation coefficient and the P-value (in brackets) as determined by the groups’ eigengene. P-values less than 0.01 are highlighted with a star. For the groups gpM1–gpM4, which exhibited high positive/negative correlations between the mRNA and protein levels, the associations were calculated according to the protein abundance. For definitions of the groups see Fig. 2.
Genes from gM5–7 were primarily expressed during the middle stage of root development (Fig. 2B). Overall, the expression of gM5 genes was higher than that of gM6 and gM7 from S1 to S3, while the expression of gM7 was higher than that of gM5 and gM6 from S5 to S7. The enriched categories included ethylene metabolism, RNA regulation transcription, and calcium signaling for gM6, protein post-translational modification for gM6 and gM7, but no enrichment was detected for gM5 (Fig. 3A). Several ethylene signaling genes (ERF1, ERF9, and ERF12) and some ethylene-responsive transcription factors (ERF72, TEM1, and RAP2.4) were found in gM6. In addition, many calcium signaling genes, including a sodium/calcium exchanger (NCL), two calcium-dependent protein kinases (CPK6 and CPK32), a calcium-binding protein (CMI1), and several calmodulin-like genes (CML3, CML5, CML23, CML36, and CML42), were also found in gM6, indicating significant roles of these genes in the middle development stage.
Genes from gM8–12 were mainly expressed at the late stages of root development (Fig. 2B). For example, gM8 and gM9 were highly expressed at S5; however, no significantly enriched categories were detected for these two groups. Genes from gM10 were highly expressed from S5 to S7, and they were enriched in amino acid (AA) degradation metabolism, cell division, development, protein synthesis, RNA processing, and RNA transcription (Fig. 3A). Several genes related to AA degradation (e.g. ASPGB1, MCC-B, and THA1) and cell division/proliferation (e.g. APC10, QQT2, and ARC3) were found in gM10. In addition, SKU5 and RHD3, which are related to root development, and RIN1, CLV2, and FCA, which are related to meristem development, were also included in this group. Genes from gM11 and gM12 had highest expression at S7. gM11 also showed relatively high expression at S5 and S6, and these genes were enriched in protein synthesis, RNA processing, and RNA regulation of transcription (Fig. 3A). No enrichment was observed for gM12.
Overall, these results revealed diverse stage-specific functions that corresponded to stage-specific expression patterns of genes exclusively regulated at the mRNA level. For example, cell wall and hormone metabolism were associated with the early stage of tuberous root development, RNA regulation transcription, ethylene signaling, and calcium signaling were associated with the middle stage, and AA metabolism, protein synthesis, and cell division were associated with the late stage.
DEGs/DEPs identified at both the transcriptomic and proteomic levels
Four groups of genes (gpM1–gpM4) were identified that were differentially expressed at both the mRNA and protein levels (Fig. 2C, Supplementary Table S6). Overall, the genes in gpM1 (Pearson’s R=0.92, P=0.003) and gpM3 (R=0.93, P=0.002) showed consistent expression patterns between the transcriptomic and proteomic levels. Group gpM1 exhibited relatively low expression from S1 to S3 but high expression from S5 to S7, with a transition at S4. Genes of this group were significantly enriched in AA metabolism, gluconeogenesis, degradation of major CHO metabolism, redox, secondary metabolism of flavonoids, and TCA/org transformation (Fig. 3A, Supplementary Table S5). Four genes involved in starch degradation were included in this group, namely ISA3, SEX4, AMY1, and PHS2. A few genes involved in TCA/org transformation (CA1, MDH, and NADP-ME3), gluconeogenesis (PCK1), and redox metabolism (MDAR2, DHAR3, GSH1, and CAT2) were also included in gpM1 (Supplementary Table S6).
In contrast, genes in gpM3 were highly expressed from S1 to S3 but had low expression from S5 to S7, again with a transition at S4. This group was enriched in AA synthesis metabolism, cell wall, major CHO metabolism, glycolysis, mitochondrial electron transport/ATP synthesis, nitrogen metabolism, redox, secondary metabolism of phenylpropanoids, and 14-3-3 protein signaling (Fig. 3A). Many genes related to the sucrose (SPP, INV1, SuSy3, SuSy4, and FK) and starch (APL3, APS1, and SBE2.1) pathways were found in this group. Several genes related to glycolysis (GAPC1, PEPC1, PK, and TPI), TCA/org transformation (NADP-ME4, NAD-ME1, and mtLPD1), and redox metabolism (CB5-B, OXP1, GR1, TPX1, and TRX2) were also in this group (Supplementary Table S6). Nitrogen metabolism is closely related to carbohydrate metabolism, and four key enzymes involved in nitrogen metabolism were observed in this group, namely NIR1, GLN1;1, GLT1, and GDH1. gpM3 also included several genes involved in mitochondrial electron transport/ATP synthesis (AOX2, COX6B, and ATP5), many genes related to the cell wall (FLA8, XTH5, SKS5, PME3, and RHM1) and phenylpropanoid metabolism (PAL, C4H, C3H, COMT, 4CL, CCoAOMT, and CAD), and some genes involved in 14-3-3 protein signaling (GRF2 and GRF11) (Supplementary Table S6).
In contrast to gpM1 and gpM3, genes in gpM2 (R=–0.87, P=0.011) and gpM4 (R=–0.93, P=0.002) showed opposite expression patterns between the transcriptomic and proteomic profiles. The enriched categories in gpM2 included major CHO metabolism, redox, signaling in sugar and nutrient physiology, and transport of major intrinsic proteins, while in gpM4 they included protein folding, protein synthesis, protein targeting, and redox of thioredoxin (Fig. 3A). Specifically, gpM2 included many genes involved in starch biosynthesis (SS1, SS4, GBSS1, DBE1, SBE2.2, and SBE2.3) and degradation (DPE1 and PHS1.1–1.4), and a number of genes related to redox (MDAR6, DHAR2, GPX6, FSD2, CSD1, and TRXF2). gpM4 also included genes related to redox (GPX4, TRX1, and CCS). Taken together, these results suggested a major contribution of post-transcriptional regulation in starch and redox metabolisms during root development.
Overall, these results revealed highly coordinated changes in gene expression (with both positive and negative correlations) between the transcriptomic and proteomic levels. These changes in expression were also clearly stage-specific (e.g. S1–S3 versus S5–S7), suggesting distinct biological roles of genes at different developmental stages.
DEPs exclusively identified at the proteomic level
Two groups of genes (pM1, pM2) were identified that were differentially expressed at the protein level but not at the mRNA level (Fig. 2D, Supplementary Table S7). The expression of proteins in pM1 was highest at S1 followed by a gradual decrease towards S7, whereas an opposite trend was seen in pM2. Enriched categories in pM1 included cell vesicle transport, fatty acid (FA) synthesis and FA elongation, protein folding, protein synthesis and targeting, G-proteins signaling, and TCA/org transformation (Fig. 3A, Supplementary Table S5). A few key genes related to lipid metabolism (CAC3, SSI2, and ACBP4), G-proteins signaling (HopM1, MIN7, and SPK1), and TCA/org transformation (CSD2, FUM2, and LPD1) were also found in this group. Genes in pM2 were enriched in AA synthesis metabolism, lipid metabolism, protein targeting, and redox (Fig. 3A). Specifically, they included several genes related to lipid metabolism (FAB1, LACS9, and ACP2), redox signaling (MDAR1, APX1, GR2, ty2, and NTRC), and AA metabolism (CBL, ATP-PRT2, ASP5, and CYSD) (Supplementary Table S7).
Taken together, the integrative analysis of the transcriptomic and proteomic data revealed highly dynamic, coordinated, and stage-specific changes in mRNAs and protein abundances, indicating a complex regulation of expression (e.g. transcriptional versus post-transcriptional) during the development of the tuberous roots. The key pathways and genes identified through our RNA-seq and iTRAQ integrative analysis are summarized in Fig. 4.
Summary of pathways/metabolites identified by multi-omics integrative analysis as being significantly affected during specific stages of tuberous root development in cassava. Development is divided into early, middle, and late stages (see Fig. 1). The genes are regulated at the transcriptomic and/or proteomic levels.
Metabolomic profiling during root development
A total of 489 metabolites were identified (Supplementary Table S8), of which 378 were significantly altered during the development of the tuberous roots. The most abundant metabolites belonged to flavones (77), lipids (53), and AA derivatives (47), followed by phenolic acids (19), carbohydrates (14), and anthocyanins (9).
Six groups of metabolites (mM1–mM6) were identified on the basis of the changes in their abundance (Fig. 2E,Supplementary Table S9). Overall, metabolites belonging to the same compound class or participating in the same biochemical pathway tended to cluster in the same group, for example 63.5% of flavone-class metabolites were grouped in mM1, of 63.6% carbohydrates and 61.5% of AA derivatives were grouped in mM5, and 46.2% and 50.0% of lipid-class metabolites were grouped in mM1 and mM4, respectively (Supplementary Table S9). Similar to the results for mRNAs and proteins, the abundances of metabolites also exhibited stage-specific patterns during root development: mM1 and mM2 clearly showed higher abundance at the early stage (S1–S3), mM4 at the middle stage (S4), while mM3, mM5, and mM6 were higher at the late stage (S5–S7, Fig. 2E). This suggested a possible biochemical connection between gene/protein expression and metabolite changes. The patterns of metabolites at the different developmental stages of the tuberous roots are summarized in Fig. 4.
Integrative transcriptomic, proteomic, and metabolomic analysis
The associations between gene/protein expression and metabolite profiling were further examined (Fig. 3B).
The metabolite groups mM1 and mM2 showed relatively high abundance at the early stage of root development, and they were both significantly positively correlated with the expression of several groups at the transcriptomic and/or proteomic levels. For example, mM2 was significantly correlated with gM1 (R=0.88, P=0.01) and gM2 (R=0.93, P=0.002) at the transcriptomic level, with gpM3 (R=0.94, P=0.002) and gpM4 (R=0.90, P=0.005) at the transcriptomic and proteomic levels, and with pM1 (R=0.93, P=0.003) at the proteomic level. This suggested that the changes in metabolites were a combined effect of genes/proteins under complex regulation at multiple layers (e.g. transcriptomic and/or proteomic).
This conclusion was supported by our observations for mM5 and mM6, which were also positively correlated with multiple groups at the transcriptomic and/or proteomic levels but which exhibited relatively high metabolite abundance at the late stage of development. Notably, mM6 was significantly positively correlated with gM11 (R=0.97, P=0.0002) at the transcriptomic level, with gpM1 (R=0.88, P=0.008) and gpM2 (R=0.88, P=0.009) at the transcriptomic and proteomic levels, and with pM2 (R=0.92, P=0.004) at the proteomic level.
On the other hand, very few significant correlations were observed for mM3 and mM4, which showed relatively high metabolite abundance at the middle stage of development. Indeed, mM3 and mM4 were only significantly positively correlated with gM9 (R=0.92, P=0.003) and gM6 (R=0.93, P=0.003), respectively, at the transcriptomic level.
This integrative analysis established possible connections between metabolites and genes/proteins, allowing us to further characterize their roles in specific metabolism pathways.
Regulation of sucrose/starch pathways
Many genes involved in sucrose degradation and starch biosynthesis, including INV1, SuSy3, SuSy4, FK, APL3, APS1, and SBE2.1, were highly expressed from S1 to S3 and then transitioned at S4 to have low expression from S5 to S7, at both the transcriptomic and proteomic levels (Fig. 5, Supplementary Table S10). In contrast, four essential genes involved in starch degradation, namely SEX4, AMY1, ISA3, and PHS2, showed the opposite expression pattern (Fig. 5A, B), and three crucial metabolites involved in the sucrose/starch pathways (sucrose, F6P, and G6P) also changed accordingly (Fig. 5C). These results strongly suggested that starch biosynthesis was mainly activated at the early stage of root development and was very limited at the late stage, and that stage S4 was a crucial period for this transition.
Changes in expression of genes and abundances of proteins and metabolites related to sucrose/starch metabolism during specific stages of tuberous root development in cassava. (A) Summary of pathways of sucrose/starch metabolism. Proteins are coded by different colors: red indicates proteins that exhibited a positive correlation with mRNA during the early stage of development; blues indicates a positive correlation during the late stage of development, and green indicates a negative correlation. (B) Heatmap of expression of genes/proteins related to sucrose/starch metabolism during root development. (C) Heatmap of three key metabolites of sucrose/starch metabolism during development. For definitions of the groups see Fig. 2. Samples/groups are identified as follows: g, transcriptomic (mRNA) level; p, proteomic level; and m. metabolomic level.
Notably, of seven DEGs identified in the production of starch from ADP-Glc, SBE2.1 was the only one that showed a consistent change at the mRNA and protein levels (Fig. 5B), suggesting that it was a key gene for this process under transcriptional regulation. The other six DEGs (SS1, SS4, GBSS1, DBE1, SBE2.2, and SBE2.3), together with two genes involved in starch degradation (DPE1 and PHS1.1–1.4), exhibited negative correlations between their expression at the mRNA and protein levels, indicating a major involvement of post-transcriptional regulation in the process.
Regulation of lignin biosynthesis pathways
Lignin is one of the main components of the plant cell wall and provides mechanical strength. There are four main branches in the pathway for lignin biosynthesis (Fig. 6A). Our metabolomic analysis identified a total of nine metabolites in lignin biosynthesis that showed significant changes during the development of tuberous roots, namely cinnamic acid, coumaric acid, caffeic acid, ferulic acid, coniferylaldehyde, coniferyl alcohol, sinapic acid, sinapaldehyde, and sinapyl alcohol. These all exhibited high abundance at the early stage and then gradually decreased towards the late stage. Correspondingly, seven DEGs involved in lignin biosynthesis were also identified, namely PAL, C4H, C3H, COMT, 4CL, CCoAOMT, and CAD, most of which were highly expressed at the early stage and then had lower expression at the late stage at both the transcriptomic and proteomic levels (Fig. 6B). The changes in the abundance of the metabolites were correlated with changes in the expression of cell wall genes (Supplementary Fig. S2). It is notable that these metabolites were primarily involved in two of the four branches of lignin biosynthesis, namely those producing coniferyl alcohol and sinapyl alcohol, indicating that these are the main pathways in cassava tuberous roots. However, a few genes, such as PAL and 4CL, showed an opposite expression pattern between the mRNA and protein levels (Fig. 6B), indicating the possible involvement of post-transcriptional regulation.
Changes in expression of genes and abundances of proteins and metabolites related to lignin biosynthesis during specific stages of tuberous root development in cassava. (A) Summary of pathways of lignin biosynthesis. Heatmaps are shown where the abundance of the metabolite changed significantly between the seven sampling times during development (see Fig. 1). Genes that were identified as being differentially expressed are indicated. (B) Heatmap of lignin biosynthesis genes/proteins that were significantly affected during specific stages of tuberous root development. For definitions of the groups see Fig. 2. Samples/groups are identified as follows: g, transcriptomic (mRNA) level; p, proteomic level; and m. metabolomic level.
Thus, the results suggested that transcriptional and post-transcriptional regulation of lignin biosynthesis during the development of tuberous roots in cassava mainly acted on the coniferyl alcohol and sinapyl alcohol branches of the pathway.
Regulation of flavonoid biosynthesis pathways
Flavonoid biosynthesis is closely associated with lignin biosynthesis as phenylalanine is the common precursor of the two pathways. A total of 16 metabolites related to flavonoid biosynthesis were identified that showed significant changes during the development of tuberous roots (Fig. 7A). Similar to the changes observed in lignin biosynthesis, they mostly showed high abundance at the early stage, although some of them were also highly accumulated at S5 (e.g. luteolin and trifolin). However, our analysis only identified three flavonoid biosynthesis genes, of which F3H and F3´H were differentially expressed at the transcriptomic level, while CHI was differentially expressed at both the transcriptomic and proteomic levels (Fig. 7B), suggesting the involvement of both transcriptional and post-transcriptional regulation in flavonoid biosynthesis in the tuberous roots.
Changes in expression of genes and abundances of proteins and metabolites related to flavonoid biosynthesis during specific stages of tuberous root development in cassava. (A) Summary of pathways of flavonoid biosynthesis. Heatmaps are shown where the abundance of the metabolite changed significantly between the seven sampling times during development (see Fig. 1). Genes that were identified as being differentially expressed are indicated in red. (B) Heatmap of flavonoid biosynthesis genes/proteins that were significantly affected during specific stages of tuberous root development. F3H and F3´H were non-significant at the protein level. For definitions of the groups see Fig. 2. Samples/groups are identified as follows: g, transcriptomic (mRNA) level; p, proteomic level; and m. metabolomic level.
Regulation of amino acid metabolism
A total of 11 AA metabolites showed significant changes, of which all except lysine increased in abundance during root development (Supplementary Fig. S3). This suggested that AA accumulation mainly occurred at the late developmental stages.
This conclusion was supported by the mRNA and protein expression profiling, as AA metabolism genes were significantly enriched in the gM10, gpM1, and pM2 groups (Fig. 3A), and they showed similar trends to the corresponding abundances of AA metabolites. In gM10, ASB1 catalyses the first step of tryptophan biosynthesis while ASN3 encodes a crucial enzyme for asparagine synthesis (Supplementary Table S4). In gpM1, APG10 encodes a BBMII isomerase that is involved in histidine biosynthesis whilst SDH encodes a saccharopine dehydrogenase that is responsible for lysine degradation (Supplementary Table S6). In pM2, CBL encodes a key enzyme in the methionine biosynthetic pathway, ATP-PRT2 catalyses the first step of histidine biosynthesis, and ASS1 encodes an arginosuccinate synthase that is crucial for arginine biosynthesis (Supplementary Table S7).
Discussion
Regulation of starch biosynthesis and degradation during tuberous root development
Transcriptional and post-transcriptional regulation are important controllers of starch biosynthesis and degradation (Smith et al., 2004; MacNeill et al., 2017), but their regulatory mechanisms remain largely unclear in cassava despite a number of studies being conducted. Chen et al. (2015) demonstrated that miR394 targets SPS2F and APL2, which participate in starch biosynthesis, whilst Zhou et al. (2019) found that that post-transcriptional silencing of SBE2 results in very high amylose starch content in the storage root. Ding et al. (2019) showed that two starch biosynthesis genes, AGPase and SS, are significantly suppressed by drought at the mRNA level but not at the protein level. Many other genes involved in starch metabolism in cassava have also been identified, but information regarding regulatory mechanisms (e.g. transcriptional and post-transcriptional) is lacking as the studies were performed only at the transcriptomic or proteomic level without a multi-omics integrative analysis (Yang et al., 2011; Wang et al., 2016b). In our current study, three genes related to starch biosynthesis (APL3, APS1, and SBE2.1), five genes related to starch degradation (BAM1, SEX4, AMY1, ISA3, and PHS2), and five sucrose metabolism genes (SPP, INV1, SuSy3, SuSy4, and FK) showed consistent expression profiles between the mRNA and protein levels during tuberous root development in cassava (Fig. 5). Moreover, their expression patterns were highly coordinated with changes in several ATP/ADP and sugar transporters (NTT1, NTT2, SWEET1, and SWEET6; Supplementary Table S11) and changes in sucrose, F6P, and G6P (Fig. 5C), and this was in line with starch accumulation in the roots (Fig. 1C). These results strongly suggested that transcriptional regulation was the primary factor controlling starch metabolism during root development, which is consistent with a previous proposal that transcriptional regulation provides mid- to long-term adjustment for starch turnover (e.g. storage starch) (Stitt and Zeeman, 2012).
Several additional genes categorized as starch biosynthesis (SS1, SS4, GBSS1, DBE1, SBE2.2, and SBE2.3) and degradation (DPE1 and PHS1.1–1.4) were also identified, but they showed significant negative correlations between the mRNA and protein levels (Fig. 5B), indicating that post-transcriptional regulation also played s significant role in starch metabolism during root development. It was notable that these genes were closely located to the pathways of sugars that are carbon links between sucrose and starch metabolism (e.g. Glc and G1P), suggesting they had a potential function in the coordination of biosynthesis and degradation of transient starch, because transient starch is used for temporary carbon storage and can be rapidly switched between the sucrose and starch pathways (Pena-Valdivia and Ortega-Delgado, 1991; Skryhan et al., 2018).
Overall, our results suggested that both transcriptional and post-transcriptional regulation played important roles in starch metabolism during tuberous root development, and probably with different emphasis; for example, transcriptional regulation for storage starch and post-transcriptional regulation for transient starch.
The early stage of development is a key period for starch biosynthesis and pathways related to the cell wall
In addition to starch biosynthesis, genes related to the cell wall, lignin, redox, energy metabolism, and signaling transduction were significantly enriched at the early developmental stages (Fig. 3). To characterize the relationships among these pathways and to further examine the regulatory networks of starch biosynthesis, APL3, APS1, and SBE2.1 were selected as query genes for co-expression analysis (Fig. 8,Supplementary Table S12).
Co-expression network of genes involved in starch biosynthesis of tuberous roots in cassava. APL3, APS1, and SBE2.1 were selected as query genes and the analysis was performed on samples taken throughout the development of the roots (S1–S7, see Fig. 1). The highly co-expressed genes are mainly associated with the functional categories carbohydrate metabolism, energy metabolism, redox, cell wall, secondary metabolism, lipid, and signaling transduction. The sizes of the nodes are based on the number of connections to other nodes.
Two sucrose degradation genes, FK and SuSy4, were highly co-expressed with APL3 and SBE2.1, indicating their significant roles in starch accumulation and root tuberization, as previously reported in potato (Baroja-Fernández et al., 2009). Nitrogen metabolism is highly related to sucrose and starch metabolism and NIR1, which encodes a key enzyme in nitrate assimilation, and AQP1, which acts as an ammonium transporter (Jahn et al., 2004), were included in the co-expression network. Glycolysis is closely linked to sucrose/starch metabolism and plays important roles in storage root development (Yang et al., 2011), and in accordance with this we identified IPGAM1, PK, and TPI, which participate in glycolysis. In addition, a few genes related to mitochondrial electron transport/ATP synthesis (ATP5) and TCA/org transformation (ACLB2, NAD-ME1, and mtLPD1) were also included in the network. These results suggested that there was a considerable demand for energy during the early development of the tuberous roots, which is in accordance with significant enrichment of glycolysis and energy-production proteins during starch accumulation (Wang et al., 2016b). Five redox genes were also included in the network, and TRX2, TPX1, and TXND9 in particular were highly co-expressed with the three genes that we queried, suggesting that they might function in starch biosynthesis via redox regulation (Geigenberger, 2003; Skryhan et al., 2018).
Cell wall genes/proteins change significantly during tuberous root development in cassava (Yang et al., 2011; Wang et al., 2016b). Three cell wall genes were found in our co-expression network, of which SKS5 encodes a protein similar to that of SKU5, which is involved in directional growth of the root tip (Sedbrook et al., 2002). Lignin metabolism is coordinately regulated with cell wall development (Gray et al., 2012). In addition, both the lignin and flavonoid pathways are highly related to starch accumulation in early storage root development (Wang et al., 2016a). Here, two genes related to the lignin pathway (C3H and COMT) and two related to the flavonoid pathway (BAN and TT5) were included in the co-expression network. These lignin biosynthesis genes were significantly down-regulated during root development (Fig. 6), which indicates similar roles of the lignin pathway in root development between cassava and potato (Firon et al., 2013).
Hormones such as auxin and ethylene play important roles in cell wall expansion and cell elongation (Majda and Robert, 2018) as well as in tuberization (Vreugdenhil and Dijk, 1989). Although no auxin biosynthesis genes were found in our co-expression network, an auxin transport, PGP1, which mediates cellular efflux of IAA and interacts with PIN genes, was included, indicating a contribution of auxin efflux but not auxin biosynthesis to the process of starch metabolism (Sojikul et al., 2015). An ethylene biosynthesis gene, ACO1, which is a crucial factor in tuberous root development in cassava (Sojikul et al., 2015), was also included. Changes in lipids alter cellulose contents and cell elongation that result in changes in the lengths of root hairs and primary roots (Li et al., 2011). Two lipid genes, FQR1 and DGK7, which are related to auxin responses and root elongation, respectively (Laskowski et al., 2002; Gómez-Merino et al., 2005), were included in our co-expression network (Fig. 8).Also included were two genes encoding 14-3-3 proteins (GRF2.1 and GRF2.2), one calcium signaling gene (CRT1a), and three G-protein signaling genes (RA-5, RAB11c, and MYA2), suggesting that they might play similar roles in cassava root development to those reported in other plants (Alvarez et al., 2011; Upadhyaya et al., 2013; Wang et al., 2016b).
The middle stage of root development is a key transition from elongation to tuberization
Jasmonic acid (JA) and Ca2+ signaling are known to participate in root tuberization in potato (Raíces et al., 2003; Nam et al., 2008). Suppression of lipoxygenase (LOX), a key protein in JA biosynthesis, dramatically disrupts tuber enlargement as well as tuber yield and shape in potato (Kolomiets et al., 2001), whilst increases in tuber number and size are positively correlated with elevated expression of LOX and two Ca2+-dependent genes, CaM1 and CDPK1 (Upadhyaya et al., 2013). CDPK is known to plays an important role in the initiation of storage roots in cassava (Sojikul et al., 2010). In our study, LOX3 and a number of Ca2+-signaling genes (including CML23, CML38, CML42, CDPK6, and CDPK32) were highly expressed specifically at the S4 middle stage of root development (Fig. 2B, group gM6), indicating that JA and Ca2+ signaling might play similar roles in cassava root development and that stage S4 is a crucial period for root tuberization.
In addition to promoting root tuberization, JA has also been shown to inhibit primary root growth (Huang et al., 2017). Over-expression of JA-responsive genes reduces root length and results in a stunted growth phenotype in Arabidopsis (Ellis et al., 2002). JAZs and MYC2 are important regulatory components of the JA signaling pathway, and they interact to mediate the inhibition of primary root growth by JA (Fernández-Calvo et al., 2011). Interestingly, a few homologues of JAZ1, JAZ3, and MYC2 were found in group gM6 (Supplementary Table S4), which exhibited highest expression at S4, thus indicating their possible roles in the inhibitory effect in root elongation during this period.
Ethylene also plays an important role in root growth by inhibiting cell expansion in the root elongation zone (Van de Poel et al., 2015). Specifically, a key transcription factor in ethylene signaling, EIN3, can interact with JAZs to positively regulate the JA-dependent inhibition of primary root growth (Zhu et al., 2011). In our study, several ethylene-response factors (including ERF1, ERF4, ERF9, ERF12, and EDF1), which are genes regulated by EIN3 (Alonso and Stepanova, 2004), had highest expression at S4, consistent with the inhibition of root elongation at this stage.
Apart from JA and ethylene signaling, lipids such as lysophosphatidylcholine (LPC) and lysophosphatidyle thanolamine (LPE) also affect primary root growth by inhibiting elongation (Li et al., 2011). Our metabolomic data indicated that many types of LPC and LPE were specifically found in group mM4 and showed greatest abundance at S4 (Fig. 2E, Supplementary Table S9). Moreover, the Ca2+-dependent lipid-binding protein CLB, the mutants of which show significantly reduced root length (de Silva et al., 2011), and several lipid metabolism genes (e.g. FAD6, KCS11, pPLAIIIb, and TAG1) had highest expression at S4, which is also consistent with root elongation being restrained at this stage.
Taken together, the expression profiles of genes related to Ca2+ signaling, JA, and ethylene signaling, and the metabolite changes observed for LPC and LPE strongly suggested that the S4 middle stage of development is a key period for the root transition from elongation to tuberization, and is in accordance with the subsequent continuation of tuberous root growth by increases in width but not in length (Lowe et al., 1982).
Conclusions
This study is the first attempt to investigate tuberous root development in cassava by means of a multi-omics integrative analysis. Our results not only reveal highly dynamic, coordinated, and stage-specific profiles of mRNAs, proteins, and metabolites during root development, but also highlight the pathways/genes that function specifically at different developmental stages. We also demonstrate that both transcriptional and post-transcriptional regulation play important roles in the development of the tuberous roots. Our findings will provide a comprehensive resource for future genetic improvement of cassava. The datasets that we have generated will also serve as a useful resource for the further investigation of tuberous root development by comparison of cassava with other root crops.
Supplementary data
The following supplementary data are available at JXB online.
Fig. S1. Principal component analysis of the transcriptomic, proteomic, and metabolomic samples.
Fig. S2. Heatmaps of genes related to the cell wall and auxin.
Fig. S3. Changes in abundance of amino acid metabolites.
Table S1. Raw counts for all genes identified in this study.
Table S2. Raw data for FPKM values for all genes identified in this study.
Table S3. Raw data for abundances of the 2657 proteins identified in this study.
Table S4. Details of the DEGs identified exclusively at the transcriptomic level.
Table S5. Values for functional enrichment presented in Fig.3A.
Table S6. Details of the DEGs/DEPs identified at both the transcriptomic and proteomic levels.
Table S7. Details of the DEPs identified exclusively at the proteomic level.
Table S8. Raw data for abundances of the 489 annotated metabolites identified in this study.
Table S9. Details of the metabolites presented in Fig. 2D.
Table S10. Details of the sucrose/starch metabolism genes involved in control of transcriptional and post-transcriptional regulation.
Table S11. Expression patterns of ATP/ADP and sugar transporters at the mRNA level.
Table S12 Details of genes involved in the starch biosynthesis co-expression network presented in Fig. 8.
Acknowledgements
This research was funded by the National Key Research and Development Program of China (grant no. 2019YFD1001100), the National Natural Science Foundation of China (31771859), the Central Public-Interest Scientific Institution Basal Research Fund for Chinese Academy of Tropical Agricultural Sciences (1630052017021, 1630012019009, 1630052019023), and the Central Public-Interest Scientific Institution Basal Research Fund for Innovative Research Team Program of Chinese Academy of Tropical Agricultural Sciences (17CXTD-28). We also thank Dr Jian Hua (Cornell University) for helping with the preparation of the manuscript and for useful discussions.
Data availability
The datasets generated and analysed in this study have been submitted to the NCBI Sequence Read Archive under the accession numbers SRR10480882–SRR10480904.
References
Author notes
These authors contributed equally to this work.
Comments