Huanan No. 9, an edible Manihot esculenta Crantz variety widely cultivated in China, and the fungus Rhizophagus irregularis (Ri), strain DOAM197198 were used as the host and AMF in this study.

Selected stems of Huanan No. 9, with the same thickness and size, were divided into uniform lengths (15cm). Next, the stem pieces were soaked in a mixed solution of acetamiprid, carbendazim, and acetamiprid to remove bacteria, fungi, eggs, etc. Then they were grown in pots (diameter of 15 cm) filled with a sterilized (high temperature of 121°C and high-pressure steam for 1 h) substrate at 2/3 height of the pot ( Figure S1 ). The substrate was composed of a mixture of river sand and vermiculite at a ratio of 4:1. Prior to sown the pieces of stem in the substrate, about 800 spores (saturated in 1 g of distilled water) were spread on the surface of the substrate of the treatment group (to induce the AM symbiosis) and covered with an amount of the sterilized substrate until the stems could stand upright in the pot. An equal amount of sterilized substrate was added to the control group. In total, three groups, including two treatment groups (Ri, inoculated with the AMF) and one control group (CK, without inoculation), were formed. Each group was composed of 20 replicates. The Ri groups were Ri-3w (allowed to grow up to three weeks after inoculation) and Ri-6w (allowed to grow up to six weeks after inoculation). The pots were maintained in a greenhouse at 28°C and 16 h photoperiod (16 h light and 8 h darkness). All the pots (CK and Ri) were irrigated every week with 100 ml of low-phosphorus Hoagland nutrient solution (phosphorus content of 20 µM).

CK and Ri-3w root samples were collected three weeks after treatment, while Ri-6w roots were sampled after six weeks. The entire root system was cut off with scissors in clean Petri dishes and gently rinsed with tap water to remove the substrate. For each treatment, two groups of samples in three replications were prepared, one for microscopic observations and the other was immediately frozen-dried in liquid nitrogen, followed by storing at -80°C for metabolomics analysis and transcriptome sequencing. From our microscopic observations, the root starts to contact with the soil mycorrhizal layer and symbiosis occurs at 3 weeks, while a stable symbiosis occurs at 6 weeks. Therefore, the comparison between CK and Ri-3w is expected to show the changes of root response at the initial stage of symbiosis. The changes in root response from the initial stage of symbiosis to the stable symbiosis stage will be elucidated by comparing Ri-3w and Ri-6w samples. Samples of 3- weeks after inoculation were selected for metabolomics and transcriptomics analyses based on microscopic observations and to investigate molecular mechanisms at the initial stage of fungal hyphae colonization of cassava roots. The Ri-6w samples were taken to examine dynamic changes in the accumulation of DAMs and expression patterns of DEGs in inoculated plants.

At the sampling time, some physiological performances were evaluated, including plant height, stem thickness, root length, root weight (biomass), and leaf chlorophyll content. For the root biomass, root samples were oven dried at 60 °C for 72 hours. The total leaf chlorophyll content was measured on three fully opened leaves, with a SPAD (single-photon avalanche diode) meter on the plants. Three technical measures per leaf were conducted.

The WGA (wheat germ agglutinin) fluorescent dye staining method was used for the mycorrhizal staining directly after the samples’ collection. The root samples were introduced in a clean 50 ml centrifuge tube containing 10% KOH solution. Then they were fixed in FAA (formaldehyde alcohol acetic acid) solution for 24 h, followed by washing with distilled water. Next, root samples were placed in a water bath at 90°C for 5 minutes and then soaked again in a 10% KOH solution. After removing the KOH by rinsing three to five times with purified water, root samples were subsequently immersed in a 2% HCl solution, gently stretched using tweezers, left at room temperature for 15-30 minutes, and then washed four to five times in distilled water. Following, root samples were soaked in PBS buffer for 30 minutes; gently stretched with tweezers; and then washed to discard the PBS buffer. Finally, the root samples were placed on a clean 2 ml centrifuge tube, and a fresh mixture solution of PBS buffer and WGA488 (V:V = 1000:1) was added. The tubes were wrapped in tin foil paper and kept overnight in the refrigerator at 4°C. The following day, the symbiotic cells were observed using scanning electron microscopy (SEM, S-3000N, Hitachi Co., Ltd., Matsuda, Japan).

The widely targeted metabolomics profiling of root samples was carried out at MetWare Co., Ltd. (Wuhan, China). Briefly, 100 mg of each root sample (beforehand vacuum freeze-dried) was ground into a fine powder and dissolved in 1 mL of methanol (70%), vortex-mixed and extracted overnight at 4°C. Thereafter, extracts were centrifugated at 12,000 g for 10 minutes, and supernatants were collected and filtered with a 0.22 μm microporous membrane. All samples’ final extract was stored in a vial for subsequent metabolomics analysis. All sample extracts were mixed equally to form the quality control (QC) samples. The data acquisition system was a UPLC-MS/MS (ultra-performance liquid chromatography-mass spectrometry) platform composed of a UPLC (SHIMADZU Nexera X2, www.shimadzu.com.cn/) and a tandem mass spectrometry (MS/MS) (Applied Biosystems 4500 QTRAP, http://sciex.com/). The qualitative identification and quantification of metabolites; and multivariate analyses, including orthogonal partial least squares discriminant analysis (OPLS-DA), principal component analysis (PCA), hierarchical clustering analysis (HCA), differentially accumulated metabolites (DAMs) analysis, and functional annotation and enrichment analyses were conducted as per reported methods (Chen et al., 2013; Dossou et al., 2021; Dossou et al., 2022). Significant DAMs were filtered out at the thresholds of VIP (variable importance in projection) ≥ 1 and t-test p < 0.05.

The total RNA extraction from cassava root samples, subsequent integrity and quality checking, sequencing on Illumina Hiseq platform, construction of cDNA library, and de novo assembly were carried out as per recently described methods (Sakamoto et al., 2019; Wang et al., 2020). Then the clean reads were mapped into the cassava reference genome (version 6.1, https://phytozome-next.jgi.doe.gov/info/Mesculenta_v6_1) (Bredeson et al., 2016), with the TopHat2 software (Kim et al., 2013). The parameter of no more than one mismatch was accepted in the alignment (Kim et al., 2013).

The expression level of each gene was normalized to the number of FPKM (Fragments Per Kilobase of transcript per Million reads) using the Cufflinks 2.0 software (Trapnell et al., 2012). The DEGs were identified using the DESeq2 software at thresholds of |log2Fold Change| ≥ 1 and p-value < 0.5 (Love et al., 2014). The functional annotations of DEGs were conducted via GO (gene ontology) and KEGG (Kyoto encyclopedia for genes and genomes) pathway enrichment analyses. The Blast2GO and KOBAS2.0 programs were used for the GO and KEGG analyses, respectively (Kanehisa and Goto, 2000; Conesa et al., 2005).

The establishment of mutualistic interactions between plant roots and beneficial microorganisms is governed by a common pathway (the SYM signaling pathway) which has been well characterized in Medicago truncatula (Bonfante and Genre, 2010; Genre and Russo, 2016). To identify SYM signaling pathway-related genes in cassava, M. truncatula SYM pathway genes were downloaded from the phytozome website and used for blast analysis against the cassava genome (Bredeson et al., 2016). Further, the DEGs related to this pathway were screened out as potential candidate genes for future studies.

To validate the transcriptome data, twelve genes with varying expression patterns were randomly selected for qRT-PCR analysis. The analysis was performed on LightCycler480 (Roche, Switzerland) real-time PCR system, with ChamQ™ SYBR1 qPCR Master Mix (Vazyme Biotech, Nanjing, China) (Song et al., 2021). The EF-1α (elongation factor 1-α) gene was used as the internal control. The specific primers of each gene are listed in Table S1 . The expression level of each gene was computed via the 2¯^ΔΔCT method (Livak and Schmittgen, 2001).

Statistical analyses and graphing were conducted using GraphPad Prism v9.0.0121 (GraphPad 159 Software Inc., La Jolla, CA, USA). The data are presented as the mean ± SD, and statistical differences were achieved by t-test at P < 0.05. Heatmaps of gene expression patterns were carried out in TBtools software (Chen et al., 2020).

To confirm the establishment of AM symbiosis between the inoculated cassava roots and the fungi, we conducted microscopic observations on root samples at three and six weeks (3w and 6w) after inoculation. The results showed that symbiotic cells were already established at 3w, and proliferating fungal hyphae could be observed ( Figures 1A–D ). After 6w of inoculation, the fungal colonized a large surface of root cells, and hyphopodium could be observed ( Figure 1D ). As shown in Figures S1A, S1B , the AM symbiosis improved the physiological performances of inoculated cassava plants (Ri). Compared to the control plants (CK), the plant height and root length of Ri plants were significantly higher ( Figures 1E, G ). The stem thickness and root biomass (root dry weight) of Ri plants were slightly increased compared to CK ( Figures 1F, H ); however, the differences were not statistically significant. There was no difference in the chlorophyll contents of CK and Ri plants ( Figure 1I ).

To examine the influence of AMF inoculation on cassava plants’ metabolism, we carried out metabolomics analysis of CK (3w) and Ri (3w and 6w) root samples. In total, 1016 metabolites, including 529 and 487 at the positive and negative ions, respectively, were identified and classified ( Table S2 ). The main classes of metabolites identified were flavonoids (18.21%), followed by phenolic acids (16.83%), lipids (16.34%), amino acids and derivatives (10.43%), organic acids (8.07%), saccharides and alcohols (6.69%), alkaloids (6.50%), and nucleotides and derivatives (6.40%) ( Figure S2 ).

Hierarchical cluster analysis (HCA) and Principal component analysis (PCA) allow for exploring the variability of metabolites among samples of the same and different groups. As presented in Figure 2 , the CK, Ri-3w, and Ri-6w samples were clustered separately on the PCA and HCA plots, indicating that their metabolite profiles were very different. Some metabolites were specifically induced in the roots of Ri plants ( Figure 2B ). To confirm the observed metabolite variation, we performed an OPLS-DA analysis. The results were supportive ( Figure S2 ). We obtained strong goodness of fit (R2X > 0.643, R2Y> 0.997) and high predictability (Q2 > 0.777) of the models ( Figure S3 ).

In order to identify major metabolites involved in AM symbiosis in cassava, we performed DAMs analysis. Significant DAMs in pairwise comparison between groups were detected at thresholds of VIP ≥ 1 and p-value < 0.05. The volcano plots are shown in Figure S4 . The results indicate that there were 25 (22 up-regulated) significantly differential metabolites in pairwise comparison between CK and Ri-3w ( Figure 3A ; Table S3 ). Meanwhile, 73 DAMs, including 40 up-regulated and 33 down-regulated, were identified between CK and Ri-6w ( Figure 3A ; Table S4 ). Venn diagram analysis showed that only eight DAMs overlapped between the two pairwise comparisons ( Figure 3B ), indicating a dynamic metabolome change in AM-inoculated cassava roots. Heatmap clustering analysis revealed that 5-aminolevulinic acid, L-glutamic acid, and lysoPC 18:2 were the most significantly induced metabolites in inoculated cassava ( Figure 3C ). The classification of the DAMs between CK and Ri-6w is shown in Figure 3D . It was noteworthy that 14 (14/18) flavonoids were down-regulated, while 7 (7/9) and 9 (9/10) nucleotides and amino acids, respectively, were up-regulated at 6w ( Figure 3D ).

To explore the molecular mechanisms involved in AM symbiosis in cassava, we carried out functional annotation and enrichment analyses of DAMs. The results showed that the DAMs between CK and Ri-3w were primarily involved in secondary metabolites biosynthesis, microbial metabolism in diverse environments, ABC transporters, aminoacyl-tRNA biosynthesis, biosynthesis of antibiotics, and amino acids metabolism ( Figure 3E ; Table S5 ). In addition to these pathways, the DAMs between CK and Ri-6w were mainly assigned to carbon and sugar metabolisms ( Figure S5 ).

To gain insight into the molecular mechanisms involved in AM symbiosis in cassava, we conducted transcriptome analyses of CK, Ri-3w, and Ri-6w root samples. The summary of the transcriptome sequencing is presented in Table 1 . The raw data and clean data ranged from 5.78 to 7.17 Gb and 5.71 to 7.07 Gb, respectively ( Table 1 ). The Q20 varied from 97.45 to 97.76%, while the Q30 ranged from 92.95 to 93.68%. The clean reads were further mapped into the cassava reference genome, and the unique and total mapped reads varied from 74.18 to 90.99% and 75.71 to 92.86%, respectively ( Table 1 ). These results show the high quality of the RNA-seq data. To examine the influence of AMF inoculation on cassava transcriptome, we conducted correlation and principal component analysis (PCA) of samples ( Figure S6 ). The results revealed that the transcriptomes of CK and Ri-3w were not too much different. Samples of Ri-6w clustered separately and could be discriminated by PC1 (72.2%), suggesting a dynamic transcriptome change in inoculated cassava along with root development ( Figure S6B ).

To unveil key differentially regulated pathways, we first screened out DEGs in pairwise comparison between groups. The volcano plots of the DEGs in the comparison of CK against Ri-3w and Ri-6w are shown in Figure S7 . We detected 693 DEGs, including 659 up-regulated and 34 down-regulated between CK and Ri-3w ( Figure 4A ; Table S6 ). Between CK and Ri-6w, there were 6,481 DEGs, of which 2,637 and 3,844 were up- and down-regulated, respectively ( Figure 4A ; Table S7 ). Of them, 88 DEGs were commonly identified in the pairwise comparison between all groups ( Figure 4B ).

We selected the DEGs in the pairwise comparison between CK and the other two groups for functional annotation and enrichment analyses. The most enriched GO terms that involve DEGs between CK and Ri-3w were related to membrane and plastid ( Figure 4C ). Meanwhile, DEGs between CK and Ri-6w were mainly associated with membranes, cell walls and junctions, and extracellular regions ( Figure S8A ). KEGG annotation and enrichment analysis of DEGs between CK and Ri-3w indicated that they are primarily involved in secondary metabolites biosynthesis, carbon metabolism, steroid and diterpenoid biosynthesis, phenylpropanoid biosynthesis, pyruvate metabolism, glycolysis/gluconeogenesis, fatty acid metabolism, and nitrogen metabolism ( Figure 4D ). The DEGs between CK and Ri-6w were mainly assigned to the biosynthesis of secondary metabolites, plant hormone signal transduction, phenylpropanoid biosynthesis, MAPK signaling pathway, plant-pathogen interaction, and metabolism of diverse sugars ( Figure S8B ). The lists of all identified pathways are shown in Tables S8 and S9 .

As we above reported, the 6w samples were taken in order to examine the dynamic changes in the expression levels of DEGs identified at 3w. Therefore, we mainly focused on DEGs between CK and Ri-3w for identifying key candidate genes underlying AM symbiosis in cassava.

AM symbiosis promotes phosphorus and nitrogen uptake and metabolisms. We identified six inorganic phosphate transporters (PHT) and 20 nitrogen transporters related DEGs ( Figure 5A ; Table S10A ). Of the PHT-related DEGs, Manes_15G190400 and Manes_03G164700 were the most significantly induced ( Figure 5A ). The identified nitrogen transporters included 5 and 15 ammonium (AMT) and nitrate (NRT) transporters, respectively. Both five AMTs were up-regulated, with Manes_05G082500 and Manes_16G119600 being the most highly induced ( Figure 5A ). Eleven of the NRTs were up-regulated, among which Manes_16G113300, Manes_01G191900, Manes_17G061600, Manes_08G106700, and Manes_18G018900 were the most significantly induced ( Figure 5A ).

We also identified eight nitrogen metabolism related-DEGs, including two nitrate reductase (Manes_18G073100 and Manes_S112700) and seven glutathione S-transferase family genes ( Figure 5B ; Table S10B ). Both two nitrate reductase genes were down-regulated, while the glutathione S-transferase genes were up-regulated. Of them, Manes_01G123300 (Log2(FC) = 5.5) and Manes_12G039900 (Log2(FC) = 3.9) were the most significantly induced ( Figure 5B ; Table S10B ).

AM symbiosis influence sugar metabolism and phenylpropanoid biosynthesis. We identified 40 sugar metabolism and seven phenylpropanoid pathway-related DEGs ( Figure 6 ; Tables S10C, D ). It is noteworthily that all these DEGs were up-regulated in Ri at 3w ( Figures 6A, B ). Of the sugar metabolism-related DEGs, one beta-glucosidase (Manes_S042900) and six UDP-glycosyltransferases (Manes_S049100, Manes_01G018000, Manes_01G018100, Manes_01G017900, Manes_07G123300, and Manes_02G125300) were more than 8.5 folds significantly induced. The gene Manes_S091700 encoding a 4-coumarate-CoA ligase-like was the most significantly induced (Log2(FC) = 7.72) in the phenylpropanoid pathway ( Figure 6B ).

Besides, 49 (four down-regulated) diverse other transporter family-related DEGs were also filtered out ( Table S10E ), and their expression patterns are shown in Figure S9 . Of them, three ABC transporters (Manes_11G015600, Manes_13G104000, and Manes_11G148900), two acyl carrier proteins (Manes_03G055600 and Manes_16G087900), one oligopeptide transporter (Manes_03G142900), one ascorbate-specific transmembrane electron transporter (Manes_17G086900), one transmembrane protein (Manes_15G046500), four wall-associated receptors (Manes_02G171700, Manes_11G070400, Manes_11G070300, and Manes_11G070200), and two lipid transporters (Manes_14G109400 and Manes_S095100) were more than 6.8 folds up-regulated in Ri at 3w.

TFs and phytohormones play essential roles in AM symbiosis. Therefore, we screened out their respective DEGs and examined their expression profiles ( Figure 7 ; Table S10F ). The main differentially regulated TFs were DELLA (9 DEGs), EP2/ERF (13 DEGs), and MYBs (5 DEGs) ( Figure 7A ; Table S10F ). As shown in Figure 7A , the genes Manes_17G055300 (MYB), Manes_05G040000 and Manes_13G048400 (EP2), Manes_13G041400 (ERF), Manes_11G083500 and Manes_14G054900 (DELLA), were the most induced TFs by the AMF inoculation.

We identified 18 phytohormone-related DEGs, including 2, 7, 2, 5, and 2 ABA (abscisic acid), auxins (AUXs), cytokinin (CYTs), gibberellic acid (GAs), and strigolactones (STs), respectively ( Figure 7B ; Table S10G ). Most of them were highly induced, and only Manes_05G135800 (CYT) was down-regulated in Ri at 3w ( Figure 7B ).

To unveil key symbiosis (SYM) signaling pathway-related candidate genes in cassava, we performed a Blast search against Medicago truncatula SYM pathway genes. In total, we identified 41 genes homologous to M. truncatula SYM pathway genes, with a similarity ranging from 37.42 to 81.99% ( Table S11 ). Thereafter, we examined their expression profiles and found that 28 of them were differentially regulated at least at one-time point. Accordingly, they were selected as potential candidate genes for further study on the SYM signaling pathway in cassava ( Table 2 ). Among them, MeSTR2, MeRAM1, MeKinG2, MeEXO70, MeABCB12, MeRAD1, MeHYP4, MeKinG1, and MeCYTB561 were more than 7 folds significantly induced in Ri at 3w.

To confirm the reliability of the transcriptome data, twelve random DEGs were selected for quantitative real-time PCR analysis. The results showed that the expression patterns of these genes through the RNA-seq and qRT-PCR were consistent (R² = 0.80; Figure 8 ), supporting the high confidence level of our findings.