The 6-month-old “Wancha No.4” (Wang et al., 2018) cutting seedlings with about 5 cm height were used as tea plant (Camellia sinensis L.) materials, which were inoculated with AMF in pot culture, and AMF isolate was Rhizophagus intraradices BGC JX04B provided by Beijing Academy of Agriculture and Forestry Sciences (Chen et al., 2017a). R. intraradices is widely used in the study of the interaction between plants and AMF, which can colonize the root of most mycorrhizal plants well. The inocula, obtained by propagating the isolate with clover (Trifolium repense L.) as hosts for 3 months in the greenhouse, were ca. 40 spores per gram in spore density. The growth substrate was the mixture of autoclaved (121°C, 2 h) soils and peat (1:1, v:v). The soils were collected from the tea plant garden of Tea Research Institute, Anhui Academy of Agricultural Sciences, Huangshan, China (29°41′18″E, 118°15′33″N), and the soil chemical properties were determined as follows: pH 5.13, organic matter content 5.17 g⋅kg^–1, available nitrogen (N) 31.40 mg⋅kg^–1, available phosphorus (P) 2.52 mg⋅kg^–1, and available potassium (K) 76.20 mg⋅kg^–1. The substrate was additionally applied with 200 mg⋅kg^–1 N ((NH₄)₂SO₄), 50 mg⋅kg^–1 P (KH₂PO₄), and 100 mg⋅kg^–1 K (KNO₃).

A pot experiment with two treatments [non-mycorrhizal (C) and mycorrhizal (T)] was set up to reveal the effect of AMF on RSA of tea plant cutting seedlings. Each pot with one tea plant cutting seedling represented one biological replicate, and three replicates, namely three pots, were set up in each treatment. Each pot was added with 950 g substrate and 50 g AMF inocula (T treatment) or sterilized inocula (C treatment). To avoid the potential influence of other microbes, we also added 5 mL filtrates (25-μm filter) of AM fungal inoculum to the C treatment. All mycorrhizal and non-mycorrhizal seedlings grew in the greenhouse under natural light conditions with 22∼30°C and 70∼80% relative humidity. Plants were harvested at 6 months after transplanting.

After total seedlings were gently washed out from the soil, shoots and roots were separated and the fresh weights were recorded, respectively. The ScanMaker i800 Plus scanner (Hangzhou WSeen Detection Technology Co., Ltd., China) was used to obtain root system pictures with high resolution, followed by quantifying the indexes of RSA by LA-S (Leaf and Root Analysis System) plant image analyzer system (Hangzhou WSeen Detection Technology Co., Ltd., China). The scanning process must be accomplished quickly to avoid the interior change in roots. Then, the root tips (about 1.5 cm length) were randomly cut down (about 0.2 g) and preserved at −80°C for RNA-seq and quantitative real-time PCR (qRT-PCR) analysis. The rest of the root was sectioned into fragments with about 1 cm length and 50 fragments were randomly picked up for the measurement of mycorrhizal colonization. Other fragments were used for carbohydrate determination. Mycorrhizal staining was performed according to Phillips and Hayman (1970). Specifically, the root segments were incubated with 10% KOH (w/v) in a water bath at 60°C for 60 min, then rinsed with tap water, bleached with alkaline hydrogen peroxide (10% H₂O₂ + NH₄OH) for 15 min, acidified in 2% HCl solution for 5 min, and stained with 0.05% Trypan blue in lactoglycerol (lactic acid:glycerol: water, v:v:v = 1:1:1) at 60°C for 30 min. After decolorizing in lactic acid: glycerol (v:v = 1:1) for 24 h, each root segment was placed onto a slide and estimated under an optical microscope (OLYMPUS BX51), and then the mycorrhizal colonization was quantified according to Biermann and Linderman (1981).

After being digested with HNO₃, the shoot P contents were measured by a spectrophotometer according to Arshad et al. (2016). The carbohydrate contents (sucrose, reducing sugar, and soluble sugar) in fresh roots were colorimetrically determined according to Xue (1985) and Gao (2006) with some modifications. Briefly, 0.2 g fresh roots were homogenized with 4 mL 80% ethanol, incubated for 30 min at 90°C, and then centrifuged at 8,000 rpm for 10 min. The residues were extracted two times as described above, and all supernatants were combined for subsequent analysis. Total soluble sugar was assayed with a mixture of 40 μL supernatant, 160 μL deionized water, and 3 mL 0.2% anthrone in 95% sulfuric acid at 90°C for 15 min followed by measurement of absorbance at 620 nm. The sucrose content was determined by mixing 40 μL supernatant, 160 μL deionized water, and 100 μL 30% KOH at 90°C for 10 min, then adding 3 mL anthranone and bathing at 45°C for 30 min. After cooling, the OD₆₂₀ (optical density) was measured. Reducing sugar was assayed with a mixture of 0.3 mL supernatant, 1.7 mL deionized water, and 1.7 mL dinitrosalicylic (DNS) acid at 90°C for 5 min followed by measurement of OD₅₂₀. The carbohydrate contents were quantified according to the linear standard curves constructed with respective sugars.

The total RNA in the root tips from C and T treatments (three independent biological replicates in each treatment) was isolated with an RNA extraction kit (Accurate Biotechnology Co., Ltd., Hunan, China) according to the manufacturer’s protocol. The RNA integrity was checked with Agilent 2100 Bioanalyzer and Agilent RNA 6000 Nano Kit, and RNA concentration was determined using Denovix DS-11 spectrophotometer (Denovix Inc., United States). One aliquot was sent to Tsingke Biotechnology Co., Ltd. (Beijing, China) for RNA-seq by using the Illumina Hiseq™ 4,000 platform and 150 bp paired-end reads were generated. After removing reads containing adapters or more than 50% N and low-quality reads from the raw reads, clean reads were obtained for all subsequent analyses. The Camellia sinensis reference genome (Wei et al., 2018), directly downloaded from the Tea Plant Information Archive (TPIA),¹ was used for the mapping of clean reads by HISAT (Kim et al., 2015). For ascertaining the differentially expressed genes (DEGs) between C and T treatments, after aligning clean reads to reference sequences by Bowtie2 (Langmead and Salzberg, 2012) and calculating the expression levels of genes by StringTie (Pertea et al., 2015), DESeq2 (Love et al., 2014) was performed with Q-value ≤ 0.05. Gene ontology (GO) enrichment analysis of DEGs was implemented by the Blast2GO software (Conesa et al., 2005). Based on Kyoto encyclopedia of genes and genomes (KEGG),² the KOBAS 2.0 software³ was used to test the statistical enrichment of DEG in KEGG pathways (Xie et al., 2011).

Another aliquot of extracted RNA was used for qRT-PCR according to the previous protocol (Chen et al., 2017a). 18s rRNA (AB120309.1) was used as a reference gene for qRT-PCR analysis in tea plants (Wei et al., 2013). All primers (Supplementary Table 1) of randomly selected genes, designed by using Primer 3 software to amplify 150–200 bp fragments, were synthesized in Tsingke Biotechnology Co., Ltd (Beijing, China). The relative expression of each target gene was calculated by using 2^–ΔΔCt method (Livak and Schmittgen, 2001).

All data were presented as mean ± SE of three replicates. The t-test was performed and the correlations between the indexes of RSAs and DEGs were calculated with bivariate correlations analysis using SPSS (Statistical Package for the Social Sciences) v.25 statistical software (SPSS Inc., Chicago, IL).

The variation partitioning analysis (Peres-Neto et al., 2006) based on redundancy analysis (RDA) was performed to quantify the contributions of different categories of DEGs to the construction of RSA in tea plant cutting seedlings. Four explanatory variables (DEGs related to metabolisms of phosphorus, sugar, plant hormone, and lipid, respectively) were conducted using “varpart” function in R vegan package (Oksanen et al., 2013), and the generated fractions (appeared as percentage) of variation were explained by different variables with or without covariable. Normally, the explanatory variable with a higher percentage contributes more greatly to the response variable than the other explanatory variable with a lower percentage (Peres-Neto et al., 2006).

To find the key genes in AM fungal regulation on RSA, the correlations between DEGs and root indexes were discovered by performing RDA and visualized in CANOCO 5 software (Braak and Smilauer, 2002). Additionally, according to previous studies (Fusconi, 2014; Chen et al., 2017a), we proposed a hypothetical model, which was specified and analyzed with Partial Least Squares Structural Equation Modeling (PLS-SEM) with the support of WarpPLS (version 6.0) software (Kock, 2017), to explore the relationships between AMF inoculation, root branching, and metabolisms of phosphorus, sugar, plant hormones, and lipid.

As shown in Table 1, mycorrhizal colonization was 26.67% in T treatment, while that was not detected in the C treatment. After AMF inoculation, plant growth was slowing. The plant height, fresh weight of shoot, and root were significantly lower than those of the C treatment, while the R/S (root/shoot) ratio was slightly lower in T treatment.

AMF inoculation significantly increased sucrose content (p < 0.01) in the root compared with the C treatment, while no significant difference was found in the contents of reducing sugar and soluble sugar (Figure 1). With respect to the phosphorus content in the aboveground part, AMF significantly decreased it (594.43 mg⋅g^–1, p < 0.05), over 100 mg⋅g^–1 less than that in the C treatment (731.00 mg⋅g^–1, Figure 1).

Overall, AMF inoculation decreased the root morphogenesis, including total root length (TRL), total root projected area (TPA), total root surface area (TSA), total root volume (TV), and total length of roots with a different diameter (Figure 2 and Tables 1, 2), but not reaching significant level. Moreover, in the mycorrhizal treatment, the average root diameter was significantly reduced, almost two times less than that in the non-mycorrhizal treatment (Table 1).

The AR and LR formation were greatly induced by AMF inoculation, except for the first and second LR (Table 3). Since the overall increase of the root biomass in the control, the second LR number, second LR/first LR ratio, total LR (TLR) number were significantly increased compared with the mycorrhizal treatment, as well as the first LR number and first LR/AR ratio but with no significant difference (Table 3). In addition to the slight increase of AR and third LR number, inoculating with AMF significantly promoted the higher classes of LR formation. Specifically, the fourth LR was only found in the mycorrhizal treatment, and ratios of third LR/second LR, fourth LR/third LR, third LR/TRL, and fourth LR/TRL were significantly elevated, indicating more formation capacity of LR in contrast to the control.

RNA-seq analysis was simultaneously performed to reveal the molecular mechanisms possibly involved in the AM Fungal regulation on the RSA of the tea plant cutting seedlings. In Table 4, overview results of RNA-seq were outlined. After trimming and applying a quality filter, clean reads were obtained and the proportion of clean Q30 bases were all over 91% (Table 4). Over 81% of these clean reads for each biological repeat mapped into the Camellia sinensis genome in order to find DEGs between the mycorrhizal and non-mycorrhizal treatment. And the RNA-seq dataset, used for analysis in this article, was deposited in National Center for Biotechnology Information (NCBI) Short Read Archive (SRA)⁴ under the BioProject of PRJNA763027.

As revealed by RNA-seq, 2047 genes were significantly expressed induced by AMF inoculation, in which 1171 genes and 876 genes were up and downregulated, respectively (Figure 3A). Most DEGs were distributed in the range of 2∼4 (933, 45.58%) and 0.25∼0.5 (682, 33.32%) fold showed in Figure 3B. The upregulated DEGs by inoculating AMF was account for 57.12% of the total number. The most 10 up and downregulated DEGs were listed in Table 6 after removing those without accurate annotations. It is noteworthy that genes related to terpenoid synthase expressed significantly higher by 38.772 or 18.398 times in the mycorrhizal treatment (Table 6). We also found that pectinesterase 2 (114257568) was more than 10 times upregulated by AMF inoculation compared with the C treatment, which might promote the entrance of AMF into the cell layers and subsequent symbiosis establishment, indirectly indicating the successful inoculation.

To validate the RNA-seq data, 12 genes were randomly selected for expression analysis via qRT-PCR. The expressions of these genes identified by qRT-PCR were as the same trend as observed in RNA-seq (Supplementary Figure 1), and there was a significantly positive correlation (R² = 0.9387, Supplementary Figure 2) between RNA-seq and qRT-PCR data, indicating the RNA-seq data were credible for exploring the mechanism of phenotypic changes of root system induced by AMF inoculation.

All DEGs were assigned to GO categories and we found 42 significant GO terms, containing 20 DEGs in biological process, 13 DEGs in a cellular component, and 9 DEGs in molecular function. GO terms of each category based on the number of DEGs were shown in Supplementary Figure 3 and Supplementary Data 3. In detail, in the category of biological process, metabolic process (697 DEGs), cellular process (568 DEGs), localization (174 DEGs), response to stimulus (166 DEGs), and biological regulation (154 DEGs) were strongly highlighted. In the cellular component category, GO terms related to cell development were enriched with more DEGs, such as 528 DEGs in the membrane, 406 DEGs in membrane part, 399 DEGs in cell, and 394 DEGs in cell part. GO terms of catalytic activity (765 DEGs), binding (544 DEGs), and transporter activity (128 DEGs) in the molecular function category were enriched with more than 100 DEGs.

Meanwhile, we also conducted the KEGG pathway analysis and the top 20 KEGG pathways, based on RichFactor, were demonstrated in Supplementary Figure 4 and Supplementary Data 3. The pathways of metabolic pathways and biosynthesis of secondary metabolites were enriched with 401 and 260 DEGs, respectively, followed by the phenylpropanoid biosynthesis pathway (84 DEGs). Additionally, a total number of 45 DEGs belonged to fatty acid metabolism.

Taking the close relationship between RSA construction and P, sugar supply and plant hormone (especially auxin) into consideration, we selected the DEGs related to metabolisms of P, sugar, and plant hormone based on the gene description, namely “Blast nr” in Supplementary Data 2. Furthermore, since lipid metabolism induced by AMF is becoming a research hotspot on account of the deficient lipid biosynthesis in mycorrhiza together with the results of KEGG enrichment analysis in this article, we also selected lipid-relative DEGs. All these DEGs of different categories were visualized in Figure 4 and more corresponding information can be found in Supplementary Data 2.

In Figure 4A, the up (22 DEGs) and downregulated (20 DEGs) DEGs of P metabolism affected by AMF inoculation were account for about a half, respectively. Four DEGs (114304683, 114283038, 114301217, and 114261192), described as purple acid phosphatase, were both upregulated by inoculating AMF. Noteworthily, compared with the C treatment, the expression of phosphate transporter PHO1 homolog 1 (Gene ID: 114261688) in tea plant root tips was nearly two times higher under the T treatment. Only 114301937 (1-deoxyxylulose 5-phosphate synthase) was expressed over 10 times in mycorrhizal treatment more than that in the C treatment, while the expressions of 114279494 (type I inositol 1,4,5-trisphosphate 5-phosphatase 2 isoform X2) and 114316147 (type I inositol 1,4,5-trisphosphate 5-phosphatase 2 isoform X2) were elevated in the C treatment, exceeding 10 times compared with the AMF inoculation. In addition, soluble inorganic pyrophosphatase (114300100) was also downregulated by AMF. With respect to sugar metabolism, the expressions of 25 out of 39 DEGs were lower in T treatment than those in C treatment (Figure 4B). Interestingly, all the four DEGs pertaining to sugar transport protein (114263353, 114280123, 114287094, and 114314873) were significantly expressed lower induced by AMF inoculation than those in the control. Same results were also found in genes related to carbohydrate esterase (114296460, 114281749, and 114273321). For DEGs belonging to plant hormone metabolism, 23 out of 39 were upregulated after AMF inoculation (Figure 4C). It is notable that all auxin-related genes (13 DEGs) were upregulated in the mycorrhizal treatment, as well as 114321176 (related to salicylic acid), while DEGs, related to cytokinin (2 DEGs), gibberellin (3 DEGs), and abscisic acid (ABA, 5 DEGs), were downregulated, except for one ABA-related gene (114293269), which is described as abscisic acid-responsive (TB2/DP1, HVA22) family protein. Moreover, there were 14 DEGs connected to ethylene with eight expressing higher in the mycorrhizal treatment. Seven out of 12 DEGs related to the ethylene responsive factor were upregulated by AMF inoculation, indicating various genes’ function ways are not in the same way. As shown in Figure 4D, a total number of 67 DEGs were assigned to lipid metabolism, wherein 41 were upregulated and 26 were downregulated by inoculating AMF. Lipase (114280426) expressed over eight times higher in the T treatment than that in the C treatment. We also found four DEGs related to lipid transfer (114312621, 114297815, 114317640, and 114303154) were upregulated by AMF inoculation.

All the results of correlation analysis between DEGs of different categories and all root indexes were placed in Supplementary Data 4. We found that all these DEGs were not significantly correlated with the first LR number and TLR/TRL ratio. Considering a large number of DEGs, we further performed the RDA to seek out key genes regulating the RSA induced by AMF inoculation. The results were showed in Figure 5 and all potential key DEGs were filled with yellow in Supplementary Data 2. In detail, initially, as showed in Figure 5A, 114308503 (type I inositol 1,4,5-trisphosphate 5-phosphatase 2 isoform X2), 114308613 (probable protein phosphatase 2C 5), and 114312586 (protein-tyrosine-phosphatase MKP1-like) had positive correlations with the formation of first LR and second LR, the numbers of second LR and TLR, TSA, TPA, TRL, TV, and average diameter (AD), while 114297081 (glycerol-3-phosphate acyltransferase 5) and 114301937 (1-deoxyxylulose 5-phosphate synthase) positively correlated with AR and higher classes of LR number (third and fourth LR), which were consistent with the correlation analysis (Supplementary Data 4). In Figure 5B, 114296460 (probable carbohydrate esterase) significantly positively correlated with the second number, TLR number, ratio of second LR/first LR, TRL, TPA, TSA, and AD (Supplementary Data 4), same to the RDA analysis. Other four DEGs (114319435, 114273482, 114314486, and 114282010) contributed to the higher level of branching and AR formation (Figure 5B), while in the DEGs related to plant hormone, these were 114311763 (auxin-induced protein PCNT115-like isoform X3), 114266460 (auxin-responsive protein SAUR32-like), and 114312929 (ethylene-responsive transcription factor), as shown in Figure 5C. 114283883 and 114323824 were involved in ethylene and ABA metabolisms, respectively, showing positive correlations with lower classes of LR and common root indexes (TPA, TRL, TV, and AD). Among the 67 lipid-related DEGs, extracted genes, namely 114302782, 114287486, 114283523, 114268197, and 114280426, were all positively correlated with third and fourth LR development and AR number (Figure 5D and Supplementary Data 4).

To quantify the overall contributions of DEGs belonging to different categories, variation partitioning analysis was performed. The results showed that the explained fractions in descending order were lipid (19.69%), plant hormone (12.69%), P (12.41%), and sugar (9.43%) when considering with other three covariables (Table 5). And the total fraction explained was 19.72%. Intriguingly, getting rid of the covariables, we found that plant hormone-related DEGs played a negative role in RSA construction, and lipid-related DEGs showed the largest explained fractions (Table 5).

We further performed PLS-SEM analysis to decipher how AMF inoculation affected the root branching of tea plant cutting seedlings through the metabolisms of phosphorus, sugar, plant hormones (primarily auxin and ethylene), and lipid. The model we constructed could be assumed to be acceptable according to Kock (2017). The value of Tenenhaus GoF, the widely accepted model fit index for PLS-based path modeling value (Henseler and Sarstedt, 2013), was 0.515, larger than 0.36 (for large effect size), indicating the acceptance of the model. And other related results of “Model fit and quality indices” were also presented in Supplementary Table 3. As shown in Figure 6 and Supplementary Table 4, except for sugar metabolism with non-significant promotion (0.311), inoculating with AMF significantly promoted metabolisms related to phosphorus (0.980), ethylene (0.706), lipid (0.688), and auxin (0.704), and directly enhanced root branching (2.500). Auxin- and phosphorus-related metabolism significantly facilitated the root branching, while lipid (−0.387) and ethylene (−1.429, p < 0.001) were negatively contributed to the root branching. In this experiment, sugar slightly promoted the root branching with no significant difference. Moreover, the indirect (1.234, p < 0.001) and total effect (3.374, p < 0.001) of AMF on root branching were both significantly positive (Supplementary Table 4). Taking together, AMF regulating root branching was mainly simultaneous through the direct pathway (or unknown ways) and indirect pathway by influencing the metabolisms of P and auxin (Figure 6).