Soil bacterial communities associated with stony soils inuence the tuber size of Tetrastigma hemsleyanum

Plants grown in stony soils have better-developed root systems and higher crop yields than those grown in non-stony soils. The roles of various physical and chemical effects of stony soils on plant growth have been published, but the roles of soil microbiota and rhizosphere microbiota have not been investigated. Methods Tetrastigma hemsleyanum plants were cultivated for two years in stony soils and in the same soil from which rock-fragments had been removed. The microbiome and the tuber transcriptome were analyzed, using multiple bioinformatics methods.

Here, we hypothesized that rock fragment-mediated changes in soil microbiota will promote the tuber growth of T. hemsleyanum through affecting the microbiota of this plant. To address these hypotheses, we analyzed the microbiota changes of T. hemsleyanum grown in stony soils and corresponding soils from which stones were removed, followed by bacterial community pro ling, taxon abundance modeling, selection pressure, and co-occurrence network analyses. In addition, the effect of the two soils on the host plant transcriptomes were also analyzed to identify any association with microbiota.

Results
Soil physicochemical parameters and tuber size of T. hemsleyanum The effects of soil type (stony and non-stony) on the physicochemical parameters of the soil and on the tuber size of T. hemsleyanum were signi cant (p < 0.05) after two years of growth (Table 1 and Fig. S1 in Additional le 1). Several physicochemical parameters of stony soils, such as electrical conductivity (EC), total carbon (TC), organic matter (OM), and AP were signi cantly higher (p < 0.05) than those of nonstony soils, regardless of whether bulk soils (BS) or root-zone soils (RZS) were studied. Table 1 Soil physiochemical properties of the elds associated with non-stony and stony soils. EC, electrical conductivity; TN: total nitrogen; TC, total carbon; OM, organic matter; NH3-N, ammonia nitrogen; NO2-N: nitrite; AP, available phosphate; AK, available potassium; AS, available sulphur. Data are mean ± SD (n = 4 in BS and n = 5 in RZS samples On the other hand, the available sulfur (AS) concentration was signi cantly lower in stony than in nonstony soils (p < 0.05) ( Table 1). The pH of RZS in stony soils was signi cantly higher (p < 0.05) than that in non-stony soils, but there was no signi cant difference between the pH of non-stony and stony soils in BS samples. In addition, there was no signi cant difference between any of the physicochemical parameters between BS and RZS of the same soil type (p > 0.05). The mean tuber size of T. hemsleyanum grown in stony soils for two years was markedly larger than that in non-stony soils (Fig. S1 in Additional le 1). For example, the average length, diameter and weight of tubers grown in stony soils were 1.7-, 1.9and 5.6-fold higher than those in non-stony soils ( Fig. S1 in Additional le 1).

Differences In Soil Microbiota Between Non-stony And Stony Soils
For all 38 soil samples, a total of 1.63 million reads was obtained after quality ltering, representing 24,700 bacterial operational taxonomic units (OTUs) at 99% sequence similarity. The bacterial alphadiversity indices, such as Chao1, observed species richness, phylogenetic diversity (PD), and Shannon diversity index, in stony soils were generally higher than those in non-stony soils, with both the BS and RZS samples showing signi cant higher values (p < 0.05) in stony soils, compared with the corresponding samples in non-stony soils.
In the soil samples from tuber surface (TS), the observed species richness and PD of stony soils were also markedly higher than the corresponding values of non-stony soils ( Fig. S2A in Additional le 1). Principal coordinate analysis (PCoA) showed that the samples from non-stony and stony soils were clearly separated along the PCoA1 axis ( Fig. S2B in Additional le 1). This distinction was further con rmed via the analysis of similarities (ANOSIM) and permutational multivariate analysis of variance (PERMANOVA) analyses, indicating that the bacterial communities of stony soils were signi cantly different from those of non-stony soils (p = 0.001) (Fig. S2B in Additional le 1). Moreover, the bacterial communities of BS, RZS, rhizosphere soil (RS), and TS in stony soils also showed signi cant differences (p < 0.05) with the corresponding communities in non-stony soils ( Table 2). To explore whether the differences in microbial structure re ected the changes in bacterial community composition, we further analyzed the differences in taxonomic identity and relative abundance of the bacterial taxa for each soil. The most abundant phyla were Proteobacteria, Acidobacteria and Actinobacteria in both non-stony and stony soils ( Fig. S3A in Additional le 1). However, the phylum Acidobacteria showed a lower relative abundance in stony than in non-stony soils, especially in the BS samples, whereas the relative abundance of the phylum Actinobacteria was higher in stony than in nonstony soils, especially in the RZS and TS samples. In addition, the relative abundances of phylum WPS-2 in BS, RZS, RS, and TS samples of stony soils were signi cantly lower (p < 0.01), and the abundances of phylum Rokubacteria in these soil samples were signi cantly higher (p < 0.01 in BS; p < 0.05 in the others), compared with those in non-stony soils ( Fig. S3B in Additional le 1).
At the family level, the Methyloligellaceae and families belonging to the Actinobacteria in BS, RZS, RS, and TS samples of stony soils were signi cantly more abundant than those of non-stony soils (Welch's ttest, p < 0.05, Storey false discovery rate (FDR)-corrected) (Fig. 1). Furthermore, the relative abundances of the Caulobacteraceae and families from WPS-2 in BS, RZS, and TS (Fig. 1A, C and D), and certain families from the phylum Gammaproteobacteria in RZS, RS, and TS of stony soils (Fig. 1B-D) were markedly lower than the corresponding values of non-stony soils (Welch's t-test, p < 0.05, Storey FDRcorrected). Moreover, in BS samples, the families from the phyla Acidobacteria and Chloro exi of stony soils had a signi cantly lower abundance than those of non-stony soils (Welch's t-test, p < 0.05, Storey FDR-corrected) (Fig. 1A).
At the OTU level, there was a total of 110 OTUs, the relative abundances of which were more than 0.5% in at least one group (Fig. S4A in Additional le 1). Among them, there were 35 OTUs in BS samples which were signi cantly different between non-stony and stony soils, and, correspondingly, there were 16, 16 and 12 OTUs in RZS, RS and TS samples, respectively, based on three screening criteria ( Fig. 2A and Additional le 2). The four samples shared ve OTUs, including four OTUs belonging to the Alphaproteobacteria (Rhizobiales_Xanthobacteraceae_OTU2030, Rhizobiales_ Methyloligellaceae_OTU6287, Elsterales_OTU5022, and Elsterales _OTU5266), and one OTU from Actinobacteria (OTU7092) (Fig. 2B).
Twenty-one OTUs with signi cant differences in at least two groups included eight upregulated OTUs and 13 downregulated OTUs found in stony soils, compared with the corresponding OTUs in non-stony soils (Fig. 2B). These discriminatory OTUs can be identi ed as key drivers of variability to distinguish the nonstony and stony soils along the rst axis, according to the non-metric multidimensional scaling analyses Co-occurrence networks of non-stony and stony soils were constructed, based on Spearman's rank correlation (Fig. 3). The results indicated that the co-occurrence network of stony soils had higher numbers of edges, average connectivity, graph density, and modularity, but had shorter average path length, compared with the network of non-stony soils (Table S1 in Additional le 1). The network of stony soils had more nodes that belonged to Alphaproteobacteria (31.0% vs 30.2%), Actinobacteria (17.8% vs 14.6%), Verrucomicrobia (5.0% vs 2.7%) and Rokubacteria (1.8% vs 0.4%), and possessed fewer nodes assigned to Acidobacteria (19.5% vs 23.1%), Chloro exi (7.92% vs 8.3%), and Gammaproteobacteria (7.1% vs 10.1%) than that of non-stony soils ( Fig. 3A-B). Netshift analysis revealed that the number of signi cant associations in the communities of stony soils was more than that of non-stony soils (277 vs 224 associations, respectively), whereas only 21 out of 522 OTU associations existed in both communities (Fig. 3C). Moreover, a total of 54 potential driver OTUs were related to the changes in microbiome composition from non-stony to stony soils (Fig. 3C). Among these OTUs, 18 out of 21 signi cantly different OTUs mentioned above were also found to play important roles in changing the networks' structure (Fig. 2B, 3C and Additional le 3).
Linking the bacterial community with soil physicochemical parameters and tuber size Pearson's correlation analysis of 21 signi cantly discriminatory OTUs and soil physicochemical parameters indicated that most physicochemical parameters, especially AS, revealed a signi cant correlation with these discriminatory OTUs, with AS having the most associations with discriminatory OTUs, followed by pH and TN (Fig. 4). To evaluate whether these discriminatory OTUs were associated with tuber size of T. hemsleyanum, we further analyzed the correlations between these discriminatory OTUs and four tuber size parameters (Fig. 4). The results indicated that the OTUs belonging to Methylologenllaceae (OTU265 and OTU6287) and Acidothermaceae_Acidothermus (OTU1962) were signi cantly positively correlated with tuber size, whereas the OTUs associated with Xanthobacteraceae (OTU2030), Burkholderiaceae _Burkholderia (OTU2662) and Chloro exi (OTU7503) had signi cantly negative correlations with tuber size (Fig. 4).

The transcriptome of tubers grown in non-stony and stony soils
To identify whether different soil types can affect the transcriptome of tubers of T. hemsleyanum, total RNAs from tubers produced in non-stony and stony soils were sequenced separately on an Illumina HiSeq™ platform. The results showed that a total of 44,817,581 clean reads were obtained, with a 50.4% GC content and a 98.7% Q20 score. The clean data were then assembled into 115,777 unigenes with an average size of 550 bp (Table S4 in Additional le 1). The assembled unigenes were annotated using BLAST, based on sequence similarity searches against nine different public databases, and a total of 70,194 (60.6%) unigenes were annotated to at least one signi cant match in all public databases (Table  S5 in Additional le 1).
Correlation analysis showed that the gene expression patterns in tubers of T. hemsleyanum grown in the same soils were similar, but there was a marked difference in pro le between non-stony and stony soils ( Fig. S6 in Additional le 1). DESeq2 analysis indicated that a total of 3,853 differentially expressed genes (DEGs), consisting of 1,145 upregulated genes and 2,708 downregulated genes, were detected in tubers of T. hemsleyanum grown in stony soils, compared with those in tubers grown in non-stony soils ( Fig. 5A and Additional le 4). To understand the biological function of these DEGs, Gene Ontology (GO) enrichment analysis was performed. The results revealed that a large number of signi cant GO terms were identi ed between non-stony and stony soils, and the top 30 upregulated terms and top 30 downregulated terms were selected for plotting (Fig. 5).
The biological processes of the top 30 upregulated terms were prominently involved in chloroplast ssion, regulation of cellular respiration, xanthine catabolism, starch biosynthesis, abscisic acid (ABA) biosynthesis, and peptide biosynthesis (Fig. 5B). The top 30 downregulated terms were closely connected with GA-and ethylene-mediated signaling pathways, response to chitin and biotic stimuli, and shoot system development (Fig. 5C). Due to the large number of DEGs between non-stony and stony soils, we used Gene Set Variation Analysis (GSVA) to calculate the abundance of GO pathways related to each soil type. A total of 33 signi cantly different (p < 0.05) GO terms, mainly associated with auxin homeostasis, photosynthesis, secondary metabolism, response to various stresses, and organic substances, were found ( Fig. S7 in Additional le 1).
The key taxa in TS potentially affected tuber size through association with various GO pathways To determine the correlation between bacterial taxa and host pathways, a network based on the Spearman's correlation analysis was constructed, using the Cytoscape software (Fig. 6). The results indicated that the key taxa in TS had a large number of associations with 25 differentially expressed GO terms. The signi cantly upregulated OTUs, belonging to the Methylologenllaceae (OTU265 and OTU6287), Beijerinckiaceae_Rhodobaastus (OTU7514) and Acidothermaceae_Acidothermus (OTU1962), showed very strong positive correlations with the pathways related to IAA and ABA biosynthesis, and photosynthesis, whereas the most downregulated OTUs were negatively associated with these pathways.

Discussion
The presence of rock fragments in soil is well-known to change the physicochemical properties of soils and subsequently the development of plant root systems, and, in this way, to promote plant growth [1,5].
In the current study, we found that soils embedded with rock fragments signi cantly increased tuber size of T. hemsleyanum, compared with soils from which the rock fragments had been removed (Fig. S1 in Additional le 1). Some researchers considered that the positive effects on plant growth in the presence of rock fragments were due, in part, to improved plant nutrition from microorganism-induced rock weathering [7,27]. The stony soils provided well spatial isolation and moisture regime, which constitute a speci c physical environment for microorganism growth and adaptation, so that the microbial community of stony soils was completely different from that of ne earth, from which rock fragments has been removed [10,28]. In the current study, the microbial diversity and structure of communities in stony soils were signi cantly different from those in non-stony soils, especially with respect to BS and RZS samples (Fig. S2 in Additional le 1). The soil physicochemical properties also differed markedly between non-stony and stony soils in terms of both RS and RZS (Table 1), and the environmental variations in AS, OM, and pH within the ve replicate plots were the most probable drivers of differential microbial community composition between non-stony and stony soils, based on the Mantel tests and DISTLM analysis (Table S2-3 in Additional le 1).
The importance of AS and pH in shaping microbial community composition has been reported in several previous studies [29,30]. Bacterial diversity was usually reduced in acidic soils, and pH also had an important effect on the composition of the soil bacterial communities [31]. Organic matter inputs generally increased the microbial diversity of soils, which, in turn, accelerated organic matter turnover [32,33]. In addition, the bacterial community structures were markedly different (p < 0.05) in the BS, RZS, RS, and TS samples between non-stony and stony soils, but there was no signi cant difference (p > 0.05) within the same soil type, based on the ANOSIM analysis (Table 2), further indicating that the soil microbiota induced by stony soils largely changed the microbiota of RS and TS in T. hemsleyanum. This was consistent with the observation that soil type was a major driver of microbiota in the root environment [34].
It has been reported that the plant rhizosphere can be rapidly colonized by bacteria from the surrounding soils. In this study, the composition of the soil microbial community changed signi cantly in response to removal of rock fragments, with these non-stony soils being richer in bacterial taxa belonging to Acidobacteria, Gammaproteobacteria, Chloro exi, WPS-2, GAL15, and Elsterales (Alphaproteobacteria) in both the BS and RZS samples ( Fig. 1 and Fig. 3A). These groups are widely distributed in ne earth, most of them showing a positive correlation with pH and a negative relationship with nutrient elements, such as soil organic carbon, TN, and available phosphorus [35,36]. The genomes of Acidobacteria encode large enzyme complements for degrading complex carbohydrates, which, in turn, impacts on the carbon distribution in the soil [37]. The bacteria from Chloro exi and WPS-2 were strongly associated with physicochemical soil properties, such as pH, AS, and TN, and Chloro exi were able to scavenge organic matter like sucrose, glucose, and N-acetyl-glucosamine under anoxic conditions [38,39].
In the current study, marked correlations between the abundance of these discriminatory taxa and soil physicochemical parameters were also observed (Fig. 4), indicating that discriminatory taxa may lead to changes in soil properties once the rock fragments had been removed, thus affecting the growth of T. hemsleyanum. By contrast, stony soils possessed more speci c bacterial taxa like Actinobacteria, Rokubacteria, Rhizobiales (Alphaproteobacteria), Myxococcales, Desulfarculaceae (Deltproteobacteria), Nitrosomonadaceae (Betaproteobacteria), and Chthoniobacteraceae (Verrucomicribia) in the BS and RZS samples (Fig. 1, Fig. 2B and Fig. 3A). It has been reported that most of these speci c taxa present in stony soils show good rock-weathering abilities, like the genera Arthrobacter, Mycobacterium, and Streptomyces in Actinobacteria [40,41], Rhizobium in Alphaproteobacteria [42], and Desulfopila in Deltproteobacteria [41]. Rock-weathering bacteria are able to release trapped mineral nutrients, like phosphorus and potassium, which are essential for plant growth and development [43]. The genomes of bacteria from the Verrucomicrobia and Rokubacteria contain some of the earliest dissimilatory sul te reductases to induce sul te/sulfate reduction, and thus play an important role in the global sulfur cycle [44]. The genome of Rokubacteria also exhibits an extensive capacity to biosynthesize diverse small organic molecules that promote mineral dissolution [37]. This may explain why stony soils were richer in nutrients than non-stony soils in both the BS and RZS samples ( Table 2). Most of these discriminatory bacteria in BS and RZS also e ciently colonized the RS and TS of stony soils as PGPR, to promote the tuber growth of T. hemsleyanum ( Fig. 2 and Fig. S5 in Additional le 1). In addition, the genomes of many bacteria belonging to phyla Actinobacteria, Rokubacteria and Verrucomicrobia contain unique gene clusters that can encode polyketide and non-ribosomal peptide synthetases [45,46], indicating that these bacteria can mediate pathogen suppression, and thus provide a healthy environment for tuber growth of T. hemsleyanum in stony soils.
Endophytic and rhizosphere microbiota, such as nitrogen-xing bacteria and PGPR, have been reported to in uence tuber number and tuber yield of tuberous plants [21,22]. The rotation of potato and legumes can signi cantly increase the abundance of anaerobic nitrogen-xing bacteria in the soil, accompanied by a marked increase in the tuber yield of potato, compared to continuous potato cropping [47]. PGPR can produce various phytohormones, like IAA, CKs, and ABA, which are necessary for the initiation and development of tubers in many tuberous plants [48]. It has been reported that phytohormones have decisive roles in tuber formation of various tuberous plant. For examples, GAs have a negative effect on the induction, growth, and maturation of tubers, while CKs have a positive in uence on them [48,49]. Appropriate IAA concentrations play an important role in promoting tuber initiation and nal tuber size, at which point the mature tubers cease growth processes and transition to dormancy in the presence of ABA [49]. Some bacterial genera from the Rhizobiales, such as Methylobacterium and Mesorhizobium, have been reported to produce IAA and CKs, and thus promote seed germination and grain yield [50,51]. The bacterial genera related to Burkholderia and Bradyrhizobium (family Xanthobacteraceae) were reported to have the ability to synthesize GA-like substances [52].
In the present study, bacterial taxa belonging to the Rhizobiales (OTU265, 3801, 6287 and 7514) and the Frankiales (OTU1962) were more abundant in stony than in non-stony soils, but the abundance of taxa belonging to the Xanthobacteraceae (OTU2030) and Burkholderia (OTU2662 and 7304) were signi cantly lower in stony than in non-stony soils, especially in the TS samples (Fig. 2). In addition, these discriminatory OTUs were also the main driver taxa between non-stony and stony soils (Fig. 3C), and showed signi cantly positive or negative correlations with tuber size parameters (Fig. 4), suggesting that a large tuber grown in a stony soil may be closely associated with these classes of nitrogen-xing and phytohormone-producing bacteria.
Plant growth and development are also strongly regulated by the host's genes. Here, the transcriptome of T. hemsleyanum tubers grown in stony soils was also markedly different from that in non-stony soils (Fig. 5). The genes signi cantly upregulated in stony soils, relative to non-stony soils, were mainly associated with the biological processes of photosynthetic activity, and certain phytohormone biosynthesis pathways, such as those of IAA, ABA, and CKs, whereas the markedly downregulated genes were more abundant in the processes involved in responding to GAs, ethylene, plant organ development, organic substances, and various biotic stresses (Fig. 5C). This was consistent with the observation that IAA, ABA, CKs and sucrose are positive modulators of tuberization in potato [48], while GA, ethylene and biotic stress suppress this process [53][54][55].
It is worth noting that most of the biological processes identi ed in the host corresponded well to the potentially bacterial functions of the discriminatory taxa in TS, and this assured us that the differences in bacterial community composition in TS were re ected in the transcriptome of the host. Further support for this conclusion was drawn from the strong correlations measured between differential GO terms and the microbial composition of TS (Fig. 6). The more positive relationship between the upregulated OTUs in stony soils and the GO terms related to auxin homeostasis and photosynthesis are due, perhaps in part, to stimulation from the available nutrients and phytohormones produced by abundant nitrogen-xing bacteria and PGPR in stony soils ( Fig. 2 and Fig. 6). It has been reported that the rhizobacteria of potato accelerated tuberization through increasing lipoxygenase gene expression [22]. The increase in photosynthesis helps to accumulate more carbohydrates, like sucrose and starch, which can promote the initiation, formation and growth of tubers [48,56]. The downregulated OTUs showed more positive correlations with responses to GAs, ethylene, or biotic stresses (Fig. 6), which may reduce the tuber growth through constraining tuberization or becoming susceptible to pathogens, due to the relative absence of more antibiotic-producing bacteria in non-stony soils (Fig. 1).
The present study has clari ed that differences in the abundance of particular bacterial taxa between non-stony and stony soils are associated with gene expression and tuber growth of T. hemsleyanum. However, the hypothetical model is still incomplete, and the bacteria and functions related to plant growth in stony soils need further exploration and research, which will be the focus of our future studies.

Conclusions
In this study, we addressed the observation that, compared with soils from which rock fragments had been removed, stony soils supported a signi cantly different soil microbiota, which was associated with larger tuber size of T. hemsleyanum. Soil microbiota in stony soils had greater bacterial diversity, cooccurrence network complexity, and greater abundance of bacterial taxa belonging to Actinobacteria, Rokubacteria, Rhizobiales, Desulfarculaceae, and Chthoniobacteraceae than the corresponding microbiota in non-stony soils. The differential soil microbiota between stony and non-stony soils may be mainly driven by soil physicochemical parameters, such as AS, OM, and pH. The discriminatory bacterial taxa of soils shaped the microbiota of RS and TS, which was strongly correlated with tuber size parameters. In addition, the potential functions of the discriminatory bacterial taxa in TS corresponded closely with gene pathways of the host, like phytohormone biosynthesis, photosynthesis, and biotic stress resistance, which are crucial for the initiation and development of tubers. These results not only help us to better understand the role of stony soils in improving plant growth but also provide a reference for increasing tuber yield through the use of microbial inocula.

Field experiment
This study was performed at a well-grown broad-leaved mixed forest, located in Dalan Town (29°47 N, 121°08 E), a division of Yuyao City, Zhejiang Province, China. The climate is subtropical with periodic monsoon rains, and an average annual precipitation of 1547 mm, with high values from April to October.
The average annual temperature is 16.2°C with an extreme maximum temperature of 41.0°C, and an extreme minimum temperature of -8.0°C. Six months before the start of this experiment, the tree density of the broad-leaved mixed forest was adjusted to about 3,000 trees per hectare with approximately 70% canopy density. The soils were composed of a loamy-clay layer with an average 70% loam and 30% different-sized rock fragments (0.2-10 cm in diameter). Ten experimental plots (0.5 m apart from each other) were set up in two parallel rows, covering an area of 6 m 2 (2 m × 3 m) per plot in this eld. The rock fragments from one row were removed to retain less than 5% rock fragments, using a hand rake after deep ploughing, and these ve replicate experimental plots were designated "non-stony soils". The rock fragments from the other row were evenly distributed in each experimental plot after ploughing to the same depth, and these plots were designated "stony soils" (Fig. S1 in Additional le 1). The six-monthsold T. hemsleyanum plants, produced from cuttings, were planted into the experiment plot, 50 cm apart, in May 2017. Only standard basic management practices, such as irrigation and weeding, were performed in the process of the experiments.

Sampling
Samples were collected in May 2019 (about 2 years after planting), and three soil-root system compartments of each experimental plot, namely root zone soil (RZS), rhizosphere soil (RS), and tuber surface soil (TS), were sampled according to the description of Shi et al. [57], with minor modi cations.
Brie y, each plant of T. hemsleyanum was taken intact from the eld, using an ethanol-sterilized shovel. Soil loosely attached to the plant roots and tubers was collected by gentle shaking, and the sample designated RZS. Then, the tubers and roots were separated, using a sterile pair of scissors, and the roots were transferred to a new 50 mL Falcon tube lled with 30 mL sterile phosphate-buffered saline (PBS) solution. Soil from the roots of each plant was collected by centrifugation at 12,000 g for 10 min following vigorous vortexing, and designated RS. One medium-sized tuber from each plant was selected for the collection of surface soil, as TS. After vigorous vortexing, the tuber was rapidly dried with sterile lter paper to measure its length, diameter and weight, and was then snap frozen in liquid nitrogen and stored at -80°C for the transcriptomic analysis. The corresponding soils from ve randomly selected plants in each replicate experimental plot were pooled together and stored at -80°C until used for microbial community analysis. A total of eight bulk soil (BS) samples, four from the stony plots and four from the non-stony plots, collected from the middle of each of the two experimental plots (Fig. S1 in Additional le 1), were also stored at -80°C for microbial community analysis.

Soil physicochemical analysis
A total of 11 physicochemical parameters of BS and RZS were assayed to evaluate the soil conditions in the non-stony and stony soils. Soil samples were air-dried at room temperature (

Plant RNA extraction and Illumina sequencing
The tubers from three non-stony and three stony experimental plots were selected to perform the gene expression analysis. Total RNA from the three biological replicates of each soil type was extracted, using RNAiso Plus (Takara Biotech, Shiga, Japan) according to the manufacturer's instruction. The Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) was used to evaluate the RNA integrity. An aliquot of 3 μg total RNA per sample was used to generate the cDNA libraries using Bioinformatics and statistical analysis

Microbiome analysis
The bacterial sequences obtained from the MiSeq runs were processed and analyzed using the Quantitative Insights Into Microbial Ecology (QIIME v1.9.1) [60]. The sequences were clustered to form operational taxonomic units (OTUs) at the 99% sequence similarity cut-off, using the pick_open_reference_otus.py script. All singletons and OTUs associated with chloroplasts or mitochondria, all unassigned and unclassi ed sequences were removed from the data before downstream analysis. To correct for the unequal sequencing depth, the OTU table was rare ed at 24,700 sequences per sample to calculate α-and β-diversity estimates, and Bray-Curtis dissimilarity between samples. Principal coordinate analysis, utilizing Bray-Curtis distances, was calculated using the PCoA() function from the R package "Ape". Non-metric multidimensional scaling (NMDS) analysis was conducted with the metaMDS() function in the R package "vegan". The analysis of similarity (ANOSIM) and permutational multivariate analysis of variance (PERMANOVA), based on the Bray-Curtis dissimilarity, were analyzed by using the Past 3.0. Differential abundance of species at the family level was determined between soils using the program Statistical Analysis of Metagenomic Pro les (STAMP), with the Welch's t tests followed by Storey false discovery rate (FDR) corrections (p < 0.05) [61]. The signi cant difference between the top 110 abundant OTUs (average relative abundance > 0.5% in at least one group) between non-stony and stony soils was measured by the t-test in the MeV software at the p < 0.05 level. The discriminatory bacterial OTUs of each group were analyzed using the linear discriminant analysis (LDA) effect size (LEfSe), based on Kruskal-Wallis sum-rank test (α = 0.05, LDA > 5.0). The contributions of differential OTUs to distinguish between non-stony and stony soils were performed by Simper analysis. The phylogenetic tree was constructed by using the neighbor-joining method with MEGA6.0 using 1,000 bootstraps. The top 600 OTUs belonging to non-stony and stony soils were retained for the co-occurrence networks analysis, based on the Spearman's rank correlation analysis. The robust correlations (|ρ| > 0.8, FDR-adjusted p-value < 0.05) between OTUs were selected and exported using Gephi 0.9.2. The "NetShift" method was used to identify potential keystone 'driver taxa' observed in bacterial co-occurrence networks between non-stony and stony soils, using the online NetShift tool (https://web.rniapps.net/netshift/), as described by Kuntal et al. [62]. The 'driver taxa' were identi ed based on the NESH score, Jaccard Index and Betweenness values, which are described in detail by Kuntal et al. [62]. Mental tests based on the Pearson's correlation were used to calculate the explanation of physicochemical variables (Euclidean distance) on the succession of the bacterial communities (Bray-Curtis distance). Distance-based multivariate analysis (DISTLM) was used to evaluate the effect of physicochemical parameters on taxonomic community structure of BS, RZS, RS, and TS, according to the description by Anderson [63]. The correlations between the 21 differential OTUs and the physicochemical parameters of soils and the parameters of tuber size of T. hemsleyanum were calculated, using the pheatmap based on Pearson's correlation analysis. SourceTracker analysis was applied to estimate the sources of the bacteria in the RZS, RS, and TS samples, as described by Knights et al. [64]. A one-way analysis of variance (ANOVA) was applied to test the differences of bacterial α-diversity indexes, bacterial taxa, and soil physicochemical parameters, using SPSS 19.0 (IBM, Armonk, NY, USA).

Transcriptome analysis
Raw reads generated by the Illumina HiSeq TM platform were initially performed to obtain clean reads by removing reads containing poly-N, adapters, and various low-quality bases. Then, the clean data were assembled into transcripts, using Trinity [65] with the setting min_kmer_cov to 2, keeping all other parameters the same. To identify the proteins that had the highest sequence similarity with the assembled transcripts, the transcripts were loaded onto public databases, including NCBI non-redundant protein (NR) and non-redundant nucleotide (NT), Swiss-Prot, Gene Ontology (GO), Clusters of Orthologous Group (COG), Kyoto Encyclopedia of Genes and Genomes (KEGG), Protein family (PFAM), Conserved domain database (CDD), and TrEMBL, and aligned by using the BlastX to retrieve their functional annotations with an E value threshold of 10 -5 . The DEGs were calculated by using the DESeq (2010) R package, based on the TPM (Transcripts Per Million) values. The p value was adjusted using a q value [66], and q value < 0.05 and |log2(foldchange)| > 2 were set as the thresholds for signi cantly different expression. The DEGs were clustered by using the pheatmap package in R. GO terms were ascribed by using the Blast2GO functional annotation, and all the DEGs were mapped to the GO database, and the number of DEGs for each term was calculated. The GO terms that were signi cantly enriched (FDRadjusted p-value < 0.05) in the biological process were visualized by using the "ggplot2" package in R. To measure the correlation between soil microbiota composition and the tuber transcriptome, the enriched GO pathways were identi ed by Gene Set Variation Analysis (GSVA) across the entire transcriptome data using the R package "GSVA" (v1.32.9), as described by Lavelle et al. [67]. Then, the GSVA output le was applied to produce the differentially expressed GO terms using the R package "limma". Finally, 3,853 DEGs were ascribed to 33 signi cant differentially expressed GO terms with the adjusted Q value < 0.05.
Given that the microbiota of TS was in close contact with tubers of T. hemsleyanum, only the key taxa of TS were used to build the network with differentially expressed GO terms. The absolute value of Spearman's correlation coe cient > 0. 6    RS, and TS samples between non-stony and stony soils. Only the OTUs with signi cant differences in at least two groups are shown. The phylogenetic tree of the 21 discriminatory OTUs, based on the 16S rRNA sequences, was performed, using the neighbor-joining method with MEGA6.0, and 1,000 bootstraps. The red taxa and blue taxa indicated that the OTUs were enriched and reduced in stony soils, respectively, compared to that in non-stony soil. Asterisks indicate a signi cant difference at *p < 0.05; **p < 0.01; and ***p < 0.001.

Figure 3
Page 26/28 Co-occurrence networks and their potential "driver taxa" between non-stony and stony soils. Cooccurrence networks of bacterial OTUs in the non-stony soil (A) and stony soil (B). A connection represents a strong (Spearman's rank correlation coe cient |ρ| > 0.8) and a signi cant (false discovery rate (FDR)-adjusted p < 0.05) correlation. The size of each node is proportional to the number of degrees.
(C) Potential "driver taxa" for changes observed in bacterial co-occurrence networks between non-stony and stony soils. Node sizes are mapped to the NESH score, which is developed to quantify the importance of the given bacterial OTUs in the association network. The nodes colored red indicated that these nodes had a positive increase in betweenness from non-stony to stony soils.

Figure 4
Pearson's correlation analysis between discriminatory OTUs, soil properties, and tuber size parameters of Tetrastigma hemsleyanum. The red taxa and blue taxa indicated that the OTUs were enriched and reduced in stony soils, compared to that in non-stony soil, respectively. The phylogenetic tree of the 21 discriminatory OTUs based on the 16S rRNA sequences was performed by using the neighbor-joining method with MEGA6.0 using 1,000 bootstraps. EC: electrical conductivity; TN: total nitrogen; TC: total carbon; OM: organic matter; NH3-N: ammonia nitrogen; NO2-N: nitrite; AP: available phosphate; AK: available potassium; AS: available sulfur; FW: fresh weight; L: length; D: diameter. Different asterisks indicate a signi cant difference at *p < 0.05; **p < 0.01; and ***p < 0.001. higher Rich factor indicates a greater degree of enrichment. Qvalue is the P value after multiple hypothesis testing and ranges between 0 and 1; the closer the Qvalue is to zero, the more signi cant is the enrichment.