Mycorrhizal symbiosis modulates the rhizosphere microbiota to promote rhizobia–legume symbiosis

Journal Pre-proof Mycorrhizal Symbiosis Modulates the Rhizosphere Microbiota to Promote RhizobiaLegume Symbiosis Xiaolin Wang, Huan Feng, Yayu Wang, Mingxing Wang, Xingguang Xie, Huizhong Chang, Like Wang, Jicheng Qu, Kai Sun, Wei He, Chunyan Wang, Chuanchao Dai, Zhaohui Chu, Changfu Tian, Nan Yu, Xuebin Zhang, Huan Liu, Ertao Wang PII: DOI: Reference:

S1674-2052(20)30433-0 https://doi.org/10.1016/j.molp.2020.12.002 MOLP 1052

To appear in: MOLECULAR PLANT Accepted Date: 7 December 2020

Please cite this article as: Wang X., Feng H., Wang Y., Wang M., Xie X., Chang H., Wang L., Qu J., Sun K., He W., Wang C., Dai C., Chu Z., Tian C., Yu N., Zhang X., Liu H., and Wang E. (2021). Mycorrhizal Symbiosis Modulates the Rhizosphere Microbiota to Promote Rhizobia-Legume Symbiosis. Mol. Plant. doi: https://doi.org/10.1016/j.molp.2020.12.002. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. All studies published in MOLECULAR PLANT are embargoed until 3PM ET of the day they are published as corrected proofs on-line. Studies cannot be publicized as accepted manuscripts or uncorrected proofs. © 2020 The Author

1

Mycorrhizal Symbiosis Modulates the Rhizosphere Microbiota to Promote Rhizobia-Legume

2

Symbiosis

3

Xiaolin Wang1,2, Huan Feng1,3, Yayu Wang4, Mingxing Wang1, Xingguang Xie1,5, Huizhong Chang2,

4

Like Wang1, Jicheng Qu6, Kai Sun5, Wei He7, Chunyan Wang3, Chuanchao Dai5, Zhaohui Chu6,

5

Changfu Tian8, Nan Yu2, Xuebin Zhang9*, Huan Liu4*, Ertao Wang1*

6

1

7

Excellence in Molecular Plant Sciences, Institute of Plant Physiology and Ecology, Shanghai

8

Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai 200032, China.

9

2

ro

of

National Key Laboratory of Plant Molecular Genetics, Chinese Academy of Sciences Center for

Shanghai Key Laboratory of Plant Molecular Sciences, College of Life Sciences, Shanghai Normal

University, Shanghai 200234, China.

11

3

College of Forestry Northwest A&F University, Yangling 712100, China.

12

4

State Key Laboratory of Agricultural Genomics, BGI-Shenzhen, Shenzhen 518083, China.

13

5

14

Technology Research Center for Industrialization of Microbial Resources, College of Life Sciences,

15

Nanjing Normal University, Nanjing 210023, China.

16

6

17

Taian, China.

18

7

Shanghai Hanyubio Co., Ltd, Shanghai 201201, China

19

8

State Key Laboratory of Agrobiotechnology, College of Biological Sciences, China Agricultural

20

University, Beijing 100193, China.

21

9

22

Province; Institute of Plant Stress Biology, State Key Laboratory of Cotton Biology, Department of

23

Biology, Henan University, Kaifeng 475001, China.

24

*

25

Lead Contact: [email protected]

lP

re

-p

10

ur

na

Jiangsu Key Laboratory for Microbes and Functional Genomics, Jiangsu Engineering and

Jo

State Key Laboratory of Crop Biology, College of Agronomy, Shandong Agricultural University,

Center for Multi-Omics Research, Collaborative Innovation Center of Crop Stress Biology, Henan

Corresponding authors: [email protected], [email protected], [email protected]

ABSTRACT

27

Plants establish symbioses with mutualistic fungi, such as arbuscular mycorrhizal (AM) fungi, and

28

bacteria, such as rhizobia, to exchange key nutrients and thrive. The plants and symbionts have

29

coevolved and represent vital components of terrestrial ecosystems. Plants employ an ancestral AM

30

signaling pathway to establish intracellular symbioses, including the legume-rhizobia symbiosis, in

31

their roots. Nevertheless, the relationship between the AM and rhizobial symbioses in native soil is

32

poorly understood. Here, we examined how these distinct symbioses affect root-associated bacterial

33

communities in Medicago truncatula, by quantitative microbiota profiling (QMP) of 16S rRNA

34

genes. We found that M. truncatula mutants that cannot establish AM or rhizobia symbiosis have an

35

altered microbial load (quantitative abundance) in rhizosphere and roots, and in particular that AM

36

symbiosis is required to assemble a normal quantitative root-associated microbiota in native soil.

37

Moreover, quantitative microbial co-abundance network analyses revealed that the AM symbiosis

38

impacts Rhizobiales-hubs among the plant microbiota and benefit the plant holobiont. Through QMP

39

of rhizobial rpoB and AM fungal SSU rRNA genes, we revealed a new layer of interaction, whereby

40

AM symbiosis promotes rhizobia accumulation in the rhizosphere of M. truncatula. We further

41

showed that AM symbiosis-conditioned microbial communities within the M. truncatula rhizosphere

42

could promote nodulation in different legume plants in native soil. Given that the AM and rhizobial

43

symbioses are critical for crop growth, our findings might inform strategies to improve agricultural

44

management. Moreover, our work sheds light on the co-evolution of these intracellular symbioses

45

during plant adaptation to native soil conditions.

Jo

ur

na

lP

re

-p

ro

of

26

46 47

Keywords ： arbuscular mycorrhizae; rhizobia; symbiosis; rhizosphere microbiota; quantitative

48

microbiota profiling; 16S; rpoB

49

50

INTRODUCTION Terrestrial plants associate with diverse microorganisms that include mutualistic, commensal and

52

parasitic microbes (Bulgarelli et al., 2013). For example, plants form mutualistic interactions with

53

arbuscular mycorrhizal (AM) fungi, one of the most important components of the soil ecosystem

54

(Miransari, 2011). AM fungi promote plant growth by facilitating the acquisition of scarce nutrients,

55

such as phosphorus, from the soil. Moreover, AM fungi influence plant diversity and connect plants

56

underground via a hyphal network that allows the movement of resources among coexisting plants

57

(van der Heijden et al., 2015).

of

51

During evolution, the genetic program for AM symbiosis was repurposed to enable other root

59

symbioses, such as that between legumes and the nitrogen-fixing rhizobia (Bonfante and Anca,

60

2009; Martin et al., 2017; Parniske, 2008). The co-evolution of microorganisms depends on

61

interactions, and in particular on cooperation and competition, among organisms. Empirical

62

studies have shown that AM and rhizobial symbiosis have synergistic effects on plant performance

63

(e.g., via the acquisition of phosphorus and nitrogen) (Afkhami and Stinchcombe, 2016; van der

64

Heijden et al., 2016; Wang et al., 2011), but, perhaps controversially, negative effects have also

65

been observed which indicate that prior inoculation with either rhizobia or AM fungi can suppress

66

the subsequent colonization by the other symbiont (Ballesteros-Almanza et al., 2010; Catford et

67

al., 2003; Hao et al., 2019; Shrihari et al., 2000). The region surrounding roots, called the

68

rhizosphere, represents the major niche for plant-microbe and microbe-microbe interactions.

69

Nevertheless, the interactions between and the co-evolution of AM fungi and rhizobia in the

70

rhizosphere of legume plants are poorly understood.

Jo

ur

na

lP

re

-p

ro

58

71

High-throughput sequencing has enabled relative abundance profiling (RMP) of microbial

72

populations, assessed by their bacterial 16S rRNA genes and fungal internal transcribed spacers

73

(ITSs). These analyses have revealed that the relative abundance (or frequency) of bacterial

74

microbiota and mycobiota in the root and rhizosphere is differentially altered in Lotus japonicus

75

mutants impaired in either the rhizobial symbiosis or mycorrhizal symbiosis (Thiergart et al., 2019;

76

Xue et al., 2019; Zgadzaj et al., 2016). Although RMP detects variation in the frequency of

77

symbiosis-associated microbiota, it does not provide information about the microbial load in a

78

sample, nor can it compare cross-kingdom microbial interactions between samples. RMP would not

79

detect links between the absolute abundance (or load) of microbes and quantitative features, such as

80

symbiosis performance. Recently, quantitative microbiota profiling (QMP) methods that can accurately examine total

82

microbial loads were used to quantitatively assess microbiota variation in environmental samples

83

(Guo et al., 2019; Tkacz et al., 2018) and we have used the QMP approach to discover novel features

84

in the assembly of bacterial microbiota within the rhizosphere (Wang et al., 2020). A quantitative

85

analysis of root-associated microbiota in native soil would illuminate the relationships between

86

rhizobial and mycorrhizal symbioses.

of

81

In this study, we used QMP to directly investigate the interactions between AM and rhizobial

88

symbioses in the rhizosphere of Medicago truncatula. Briefly, we grew wild type M. truncatula, as

89

well as mutants defective in AM symbiosis and/or rhizobial symbiosis, in native soil conditions and

90

quantified the microbial load in rhizosphere and root samples (Guo et al., 2019; Tkacz et al., 2018;

91

Wang et al., 2020). Besides absolutely quantifying the 16S rRNA genes for bacterial community

92

profiling, we also used QMP to calculate the microbial loads of rhizobial rpoB and AM fungal SSU

93

rRNA genes.

-p

re

lP

na

ur Jo

94

ro

87

95

RESULTS

96

Mutants Defective in AM Symbiosis Exhibit a Distinct Phenotype and Root Transcriptome in

97

Native Soil To survey how AM and rhizobial symbiosis affect root-associated microbiota, Medicago plants

99

were grown in native soil under controlled greenhouse conditions, and the microbiota of rhizosphere

100

soil (defined as soil particles closely attached to roots) and root samples (nodules were removed and

101

roots were ultrasonically cleaned in PBS buffer) of Medicago plants were surveyed (see Materials

102

and Methods). We examined wild type M. truncatula (Jemalong A17), as well as mutants that disrupt

103

AM symbiosis alone (M-R+: pt4-1), rhizobial symbiosis alone (M+R-: ipd3-1, nfp-1 and lyk3-1) or

104

both (M-R-: dmi2-1 and dmi3-1). We designated the mutants as M+/- or R+/- according to the mutant

105

phenotypes have been reported under simple and static laboratory conditions. LYK3, NFP and IPD3

106

play important roles in rhizobial symbiosis, but have a minor or no effect on AM symbiosis (Horvath

107

et al., 2011; Oldroyd, 2013; Zhang et al., 2015). More interestingly, a recent study indicated that

108

colonization of AM fungi in ipd3 ipd3l is sensitive to high-P levels, suggesting that heterogeneous

109

soil conditions may contribute to variability in AM colonization (Lindsay et al., 2019). PT4 is

110

essential for AM symbiosis only (Javot et al., 2007), whereas DMI2 and DMI3 are essential for both

111

AM and legume-rhizobial symbioses (Oldroyd, 2013; Parniske, 2008).

Jo

ur

na

lP

re

-p

ro

of

98

112

Initially, we examined AM symbiosis in the roots of these plants by sequencing and quantifying

113

the 18S rRNA gene for AM fungi. Since the 18S rRNA gene primers (see Materials and Methods)

114

used are specific to the Glomeromycota, only the microbial load (not the frequency) of AM fungi was

115

surveyed. We identified 28 Glomeromycota amplicon sequence variants (ASV) belonging to the 10

116

currently known AM virtual taxa (VT) (Opik et al., 2010) in our samples, with the Claroideoglomus

117

sp. VTX00193, and Glomus sps. VTX00114 and VTX00067 dominating in all VT (Supplementary

118

Table 1). The microbial load of AM VTs was lower in the roots of the M-R+ and M-R- mutants

119

compared to WT, as expected, and none of VT detected in the roots of arbuscular mycorrhizal

120

symbiosis dysfunction mutants of dmi3 and pt4 (Figure 1A and Supplementary Table 1). Intriguingly,

121

a comparable or slightly higher microbial load of AM VTs was found in the roots of M+R- mutants

compared to WT in native soil conditions (Figure 1A and Supplementary Table 1). An unconstrained

123

principal-coordinate analysis (PCoA) confirmed these findings (Supplementary Figure 1): while the

124

first principal coordinate (PCo1) displayed the compositional transition across the M+ and M- plants.

125

To verify the function of AM VTs, we examined the total phosphorus in the roots of these plants.

126

Consistent with the observed microbial load of AM VT in roots, the total phosphorus in roots was

127

highest in M+R- plants, followed by wild type A17 (Figure 1B). The phosphorus content was

128

significantly lower in the M-R+ and M-R- mutants compared with wild type and M+R- plants (Figure

129

1B).

of

122

We also performed RNA-seq profiling of roots to investigate AM-related gene expression, and

131

found that M- mutants (pt4-1, dmi2-1, dmi3-1), exhibited a distinctive transcriptome when compared

132

to M+ plants (WT, ipd3-1, nfp-1 and lyk3-1) (Figure 1C). Genes previously identified as

133

mycorrhizal-induced genes (Afkhami and Stinchcombe, 2016) were not induced in M- plants (cluster

134

5, 6 and 7; Figure 1C and Supplementary Figure 2A). By contrast, genes involved in the phosphate

135

starvation response (PSR), defense response, and response to oxidative stress (ROS) (cluster 3) were

136

more highly expressed in M- roots than M+ roots (Supplementary Figure 2B). Moreover, genes

137

required for AM symbiosis, including fatty acid biosynthetic genes (MtKAR, MtKAS II, MtFatM etc.),

138

the ammonium transporter MtAMT2, as well as the symbiotic signalling genes (MtDELLA1, MtDMI2,

139

MtDMI3, MtNSP1, MtNSP2, MtRAM1, MtRAM2, MtWRI5a, MtWRI5b, MtRAD1, MtDMI1 etc.),

140

were more highly expressed in M+ roots than in M- roots (Supplementary Figure 3).

141

Symbiosis Affects both the Frequency and Load of Microbiota in Rhizocompartments

Jo

ur

na

lP

re

-p

ro

130

142

Next, we surveyed the bacterial microbiota in the rhizosphere soils and roots of these plants by

143

analyzing 16S rRNA gene frequency and load. We generated two sequencing data sets by scaling the

144

read counts with respect to two different reference points: (a) relative to sequencing depths for each

145

sample, which yields the standard relative abundance of 16S rRNA genes (the standard RMP or

146

frequency, Figure 1D), and (b) relative to 16S spike-in reads (Supplementary Table 2) and 16S rRNA

147

gene copy numbers for each ASV (Supplementary Table 3), retrieved from the rrnDB database

148

(Stoddard et al., 2015), as we previously described (Wang et al., 2020). This approach allowed us to

149

extrapolate the quantitative abundance of bacterial cells (the 16S copy-number-weighted QMP or

150

load, Figure 1E) for each sample. The standard RMP revealed that WT and mutant plants have significantly different frequencies of

152

bacterial phyla in both the rhizosphere and root (Figure 1D). We observed a relative decrease of

153

Alphaproteobacteria and a relative increase of Gammaproteobacteria and Betaproteobacteria in the

154

rhizosphere samples from M- mutants compared to WT, indicating that dysfunction in AM symbiosis

155

affects rhizosphere occupation by Proteobacteria. In addition, we observed a relative decrease of

156

Actinobacteria and a relative increase of Alphaproteobacteria, Chloroflexi and Deltaproteobacteria

157

in root samples of M- mutants compared to those of WT, and a relative decrease of Bacteroidetes and

158

Betaproteobacteria in M+R- mutants vs WT roots, suggesting that both AM symbiosis and nodulation

159

influence root colonization by different subsets of bacterial taxa.

-p

ro

of

151

Our QMP data revealed that the bacterial loads in roots and rhizosphere soils were also

161

significantly different between WT and symbiotic mutants (Figure 1E). Compared to WT, M+R-

162

mutants showed an elevated bacterial load in the rhizosphere but a reduced bacterial load in the root

163

(Kruskal Wallis test; Rhizosphere: QAWT = 6.44 × 109, QAM+R- = 8.25 × 109, PFDR = 0.021; Root:

164

QAWT = 1.24 × 109, QAM+R- = 9.93 × 108, PFDR = 0.053). On the contrary, M-R- mutants had a

165

reduced bacterial load in the rhizosphere but elevated load in root (Kruskal Wallis test; Rhizosphere:

166

QAWT = 6.44 × 109, QAM+R- = 4.87 × 109, PFDR = 0.023; Root: QAWT = 1.2 × 109, QAM-R- = 2.1 × 109,

167

PFDR = 0.096), relative to WT. Interestingly, the pt4 mutant (M-R+) had the lowest bacterial load in

168

the root (Kruskal Wallis test, QAWT = 1.24 × 109, QApt4 = 6.59 × 108, PFDR = 0.047), which cannot be

169

attributed solely to reduced phosphorus content (compare Figure 1E and Figure 1B). The QMP data

170

also revealed that changes in the frequency of bacterial phyla observed by RMP did not always

171

reflect a concordant change in their load (Figure 1D and E). For example, Alphaproteobacteria,

172

Gammaproteobacteria and Betaproteobacteria displayed an increased frequency in the rhizosphere

173

of ipd3 relative to WT based on RMP, but their load was not significantly different based on QMP.

174

Notably, some bacterial phyla increased in load based on QMP but decreased in frequency based on

175

RMP, or vice versa. For example, Gammaproteobacteria had a reduced frequency in the rhizosphere

176

of nfp (M+R-) relative to WT based on RMP, but an increased load based on QMP. Similarly,

Jo

ur

na

lP

re

160

177

Alphaproteobacteria displayed an increased frequency in the root of pt4 relative to WT based on

178

RMP, but a decreased load based on QMP. (Figure 1D and E).

179

Defects in AM Symbiosis Reduce the Load of Rhizobiales Clades in the Rhizosphere Next, we examined the QMP data set to determine how AM symbiosis affects the bacterial

181

communities in the rhizocompartments, and compared the results to the standard RMP.

182

Permutational multivariate ANOVA based on Bray-Curtis distance matrices (adonis2) of the QMP

183

data set revealed a marked contribution of the plant genotype to the composition of rhizosphere and

184

root microbiota (Rhizosphere: R2 = 0.436, P < 0.0001; Root: R2 = 0.332, P < 0.0001). PCoA plots of

185

the same data set exposed a clear distinction between M+ plants versus M- plants (Supplementary

186

Figure 4A, B). CAP (canonical analysis of the principal coordinates) in combination with ANOSIM

187

(analysis of similarities) were used to better quantify the influence of plant genotype on the beta

188

diversity (Figure 2A, B and Supplementary Table 4). The bacterial assemblies in all M- mutants were

189

clearly separated from those that could form AM symbiosis (WT and M+R- mutants) along the first

190

principal coordinate, across both the rhizosphere [RANOSIM = 0.68, P < 0.0001; adj-R2CAP = 0.27, P <

191

0.0001] and root samples [RANOSIM = 0.36, P < 0.0001; adj-R2CAP = 0.10, P < 0.0001] (Figure 2A,B

192

and Supplementary Table 4). In contrast, defective rhizobial symbiosis alone (M+R-) did not lead to

193

consistent effects on the load of bacterial communities in either the rhizosphere [RANOSIM = 0.15, P =

194

0.126; adj-R2CAP= 0.052, P = 0.052] or root [RANOSIM = -0.062, P = 0.63; adj-R2CAP = 0.005, P =

195

0.303] compartments (Figure 2A, B and Supplementary Table 4). Interestingly, based on the standard

196

RMP of 16S rRNA genes, the nfp mutant (M+R-) assembled a different rhizosphere (RANOSIM = 0.484,

197

P = 0.007; adj-R2CAP = 0.114, P = 0.008) and root (RANOSIM = 0.3, P = 0.046; adj-R2CAP = 0.072, P =

198

0.023) microbiota compared with WT (Supplementary Table 4), consistent with observations in

199

Lotus japonicus (Zgadzaj et al., 2016). We conclude that AM symbiosis is required to maintain a

200

normal rhizosphere and root microbiota.

Jo

ur

na

lP

re

-p

ro

of

180

201

We further assayed the ASVs that consistently increased or decreased in load in roots and

202

rhizosphere soils of symbiosis-deficient mutants compared with WT. We observed a total of 85 ASVs

203

with an altered load in the rhizosphere bacterial communities of WT vs mutants (Figure 2C), and 52

ASVs with altered loads in root bacterial communities (Supplementary Figure 5) [Kruskal Wallis test,

205

PFDR < 0.05]. Notably, we found that 20 ASVs were consistently decreased in rhizosphere soils of M-

206

mutants (dmi2, dmi3, pt4), with 9 affiliated to Rhizobiales, including 3 Rhizobium ASVs, 2

207

Sinorhizobium ASVs, 2 Mesorhizobium ASVs and 1 Bradyrhizobium ASV (Figure 2C). We also

208

observed 5 ASVs that were consistently increased in rhizosphere soils of M- mutants, which were

209

also increased in rhizosphere soils of M+R- mutants (Figure 2C), indicating that these 5 ASVs are

210

regulated by both types of symbioses. Taken together, our results indicate that plant dysfunction in

211

AM symbiosis exerts a significant effect on the load of bacterial communities in the M. truncatula

212

root-soil interface, mainly reflecting the depletion of bacterial taxa from Rhizobiales clades in the

213

rhizosphere.

214

AM Symbiosis Impacts Rhizobiales-Hubs in the Bacterial Co-abundance Network

re

-p

ro

of

204

Emerging microbiota studies reveal that root associated microbiomes are structured and form

216

interconnected microbial networks, and microbial hubs may act as mediators between the plant and

217

its microbiomes (van der Heijden and Hartmann, 2016). To further characterize the AM symbiosis

218

effect on root-associated microbiomes, we analyzed the bacterial co-abundance network based on

219

significant correlations (Spearman's correlation test, PFDR < 0.05) between paired ASVs in the QMP

220

data set. Our results suggested 38.6% of the edges interacting with ASVs belonged to Rhizobiales in

221

the rhizosphere of WT and M+R- plants (ipd3-1, nfp-1 and lyk3-1) (Figure 3A). These data suggest

222

that Rhizobiales ASVs could be a microbial hub in the rhizosphere co-abundance network of plants

223

that can form AM symbiosis. Indeed, 2/6 of ASVs which were assigned as hubs based on their

224

degree, betweeness centrality and closeness centrality belonging to Rhizobiales (Supplementary

225

Figure 6A, B). In the rhizosphere co-abundance network of plants that cannot form AM symbiosis,

226

only 33.8% of the edges interacting with ASVs belonging to Rhizobiales (Figure 3B), and none of

227

the Rhizobiales ASVs were identified as hubs based on the three measurements of centrality

228

(Supplementary Figure 6C, D). Moreover, a higher mean degree centrality for Rhizobiales ASVs was

229

observed in the M+ rhizosphere co-abundance network relative to the M- rhizosphere co-abundance

230

network (Figure 3C, Kruskal Wallis test, PFDR = 0.009). When the Rhizobiales ASVs were removed

Jo

ur

na

lP

215

231

from the QMP data set, the AM symbiosis-explained variance for rhizosphere bacterial assemblies

232

declined to 18% from 27%, once again indicating an association of the Rhizobiales ASVs with AM

233

symbiosis (Supplementary Table 5). Rhizobiales ASVs were also identified as the main hubs in M+ root co-abundance network,

235

while Actinomycetales and Xanthomonadales ASVs act as hub taxa in the M- root co-abundance

236

network (Figure 3D-F and Supplementary Figure 6E-H). Intriguingly, the root bacterial

237

co-abundance network displayed higher connectivity in M- plants (Figure 3E; nodes, 300; edges,

238

4705; transitivity, 0.59) compared with M+ plants (Figure 3D; nodes, 264; edges, 3520; transitivity,

239

0.46). A highly connected microbiota network is thought to be associated with a stressful

240

environment, such as pathogen invasion (Carrion et al., 2019; Faust and Raes, 2012). The assay of

241

correlations between the load of hub ASVs in the M- root co-abundance network and the

242

transcriptome profiles of PSR-associated genes or the defense response genes showed that the hub

243

ASVs belonging to Streptomyces, Micromonospora and Lechevalieria genera which have been

244

recognized as important sources of antibiotics (Li et al., 2019a) were negatively correlated with PSR

245

associated genes like PHO1, SQD1, SQD2 and SPX3 (Supplementary Figure 7A, all ρ < −0.75, PFDR

246

< 0.001); as well as negatively correlated with pathogenesis−related genes like PR1 and PR2

247

(Supplementary Figure 7B, both ρ < −0.80, PFDR < 0.001). Thus, the highly connected co-abundance

248

network in the M- plants may indicate that root microbiota of M- plants are under stress, such as

249

limited phosphorus availability or potential pathogen invasion, which is also consistent with our

250

observation that PSR, defense response and ROS genes were increased in roots of M- plants

251

(Supplementary Figure 2B).

252

AM Symbiosis Promote the Accumulation of Rhizobia in the Rhizosphere

Jo

ur

na

lP

re

-p

ro

of

234

253

To further assay the potential impact of AM symbiosis on rhizobial taxa and to overcome

254

limitations in 16S rRNA analyses for discerning bacterial species, we specifically quantified rhizobial

255

species by sequencing their rpoB genes. Compared to the low resolution 16S rRNA gene, the rpoB

256

gene was shown to be very powerful for discriminating between closely related species (Case et al.,

257

2007; Vos et al., 2012). A total of 56 rpoB ASVs belonging to 16 documented rhizobial species of

258

leguminous plants were identified, with Rhizobium anhuiense (55.0%) and Pararhizobium giardinii

259

(14.9%) dominating the rhizosphere soil, and S. medicae (55.1%) and R. anhuiense (21.1%)

260

dominating the root compartment (Supplementary Figure 8 and Supplementary Table 6). Compared

261

to rhizosphere and root samples, only a few of rhizobial species were identified in bulk soil

262

(Supplementary Figure 8 and Supplementary Table 6). Intriguingly, we found that the load of rhizobial rpoB reads was significantly lower in rhizosphere

264

soil of M- mutants compared to WT or M+R- mutants [Kruskal Wallis test, PFDR < 0.05] (Figure 4A

265

and Supplementary Figure 9), and R. anhuiense, P. giardinii and Bradyrhizobium japonicum were

266

consistently decreased in the rhizosphere soil of M- mutants compared to WT [Kruskal Wallis test,

267

PFDR < 0.05] (Figure 4B and Supplementary Figure 8). S. medicae, which nodulates M. truncatula,

268

was retrieved only in the rhizosphere of WT (Figure 4B and Supplementary Figure 8). Unexpectedly,

269

the load of rhizobial rpoB reads in roots showed no significant consistent difference across plant

270

genotypes, even in rhizobial symbiotic mutants (Supplementary Figure 9), probably balanced out by

271

the M. truncatula non-compatible rhizobial species like R. anhuiense, P. giardinii and B. lablabi,

272

which are abundant, but showed no difference in roots among genotypes (Supplementary Figure 8).

273

However, the rpoB reads of M. truncatula compatible rhizobia (S. medicae pASV2 and pASV179, S.

274

meliloti pASV91 and pASV247) were only found in roots of WT, ipd3 and pt4 plants

275

(Supplementary Figure 8 and Supplementary Figure 10), indicating that plants that can develop

276

infection threads harbour a comparatively large load of symbiotically compatible rhizobia in their

277

roots [Kruskal Wallis test, PFDR = 0.002].

Jo

ur

na

lP

re

-p

ro

of

263

278

To test whether the accumulation of rhizobia by AM symbiosis occurs in non-legume plants as

279

well, we also quantified AM VTs and rhizobia in rhizosphere soil of three maize (Zea mays) inbred

280

lines (Li et al., 2019b) which show different levels of AM fungi colonization (Supplementary Figure

281

11A). Indeed, we found the total abundance of rhizobial rpoB reads (Supplementary Figure 11B) and

282

the abundance of P. giardinii, S. adhaerens and S. meliloti (Supplementary Figure 11C) were

283

significantly higher in rhizosphere soil of those lines with higher levels of AM colonization

284

(GEMS56 and 4F1) compared to lower AM-colonized line 975-12 [Kruskal Wallis test, PFDR < 0.05].

285

Together, our results revealed that AM symbiosis promotes the accumulation of rhizobia in the

286

rhizosphere.

287

The AM-Promotion of Microbial Communities in the Rhizosphere Promotes Nodulation in

288

Different Legumes To assay whether the AM-promotion of rhizobia accumulation in the M. truncatula rhizosphere

290

affects root nodulation in native soil, we examined the nodule numbers of the wild type and mutant

291

plants. Consistent with the published data (Parniske, 2008), lyk3, nfp, dmi2 and dmi3 mutants did not

292

form nodules, and ipd3 formed small white bumps in its roots. Intriguingly, compared to WT plants

293

(4.3 ± 1.36; mean ± SD), only 25% of pt4 mutants (0.4 ± 0.90) formed functional nodules in native

294

soil (Fisher's exact test, P = 0.0003), but 100% of pt4 mutants formed functional nodules when

295

grown in sterile vermiculite and inoculated with rhizobia (S. meliloti Rm2011, 5.4 ± 1.91; S. medicae

296

RS90, 6.24 ± 2.22; OD=0.001) (Fisher's exact test, P = 1; Figure 4C and D). Since the pt4 mutant has

297

never been reported to be affected in nodulation under simple and static laboratory conditions (Javot

298

et al., 2007), we re-examined the root nodulation of these plants as well as the ram1 mutant which

299

has been reported to be specifically affected mycorrhizal signaling and appears to have no role in

300

nodulation signaling (Gobbato et al., 2012). Consistent with the above results, all of the WT plants

301

(6.4 ± 4.72) formed functional nodules in native soil, but none of the pt4 plants were nodulated, and

302

only one functional nodule was observed in the ram1 plants (Supplementary Figure 12, 12

303

plants/genotype were used). The arbuscule phenotypes of these plants (Supplementary Figure 12)

304

were consistent with the result of the Glomeromycota 18S rRNA gene survey (Figure 1A and

305

Supplementary Table 1), none of the M- (dmi2, dmi3, pt4 and ram1) plants formed any functional

306

arbuscules in our native soil conditions. These data suggest that, in native soil, the promoted rhizobia

307

accumulation due to AM symbiosis may in turn promote nodulation of M. truncatula.

Jo

ur

na

lP

re

-p

ro

of

289

308

To compare the symbiotically compatible symbiont(s) in the native soil and root nodules, the

309

microsymbionts of 10 Medicago nodules were isolated and identified by rpoB sequencing, as

310

described earlier (Wang et al., 2018). Results showed that 8 of the 10 isolates were identified as S.

311

medicae and had rpoB sequence identical to the strain isolated from the rhizosphere soil (S. medicae

312

RS90), as well as to the rpoB reads retrieved from high-throughput sequencing (pASV2 and

313

pASV179) (Supplementary Figure 13). Only one isolate (microsymbiont7) was identified as S.

314

meliloti (the rpoB sequence had 100% identity to pASV91 and pASV247) (Supplementary Figure

315

13), the remaining one (not shown) was identified as Enterobacter sp. To assay whether AM symbiosis-promoted microbial communities of the M. truncatula

317

rhizosphere are physiologically important for leguminous plants in native soil, we analyzed the

318

nodulation of four legumes. Specifically, we inoculated Vicia faba, Pisum sativum, Phaseolus

319

vulgaris and Glycine soja with diluted rhizosphere microbial communities from M. truncatula WT

320

and mutant plants (lyk3, dmi3 and pt4). None of V. faba plants were nodulated. For the other three

321

legumes, the number of functional nodules per plant was significantly decreased when inoculated

322

with diluted rhizosphere soil from M- plants (pt4 and dmi3) compared to WT plants (Figure 4E,

323

Kruskal Wallis test; all PFDR < 0.05). In P. sativum plants, the diluted rhizosphere soil of M- plants

324

induced 70–85% fewer nodules (pt4 rhizosphere, 2.4 ± 0.42, mean ± SD; dmi3 rhizosphere, 1.3 ±

325

0.29) than WT plants (8.61 ± 1.38) (Figure 4E). Nodule numbers were similarly decreased (75-90%)

326

for G. soja plants (WT rhizosphere, 2.28 ± 0.41; pt4 rhizosphere, 0.56 ± 0.18; dmi3 rhizosphere, 0.22

327

± 0.10) (Figure 4E). In P. vulgaris plants, however, nodulation was only observed in roots inoculated

328

with diluted rhizosphere soil of WT and lyk3, whereas no nodules were observed following

329

inoculation with diluted rhizosphere soil of pt4 and dmi3 (Figure 4E). Fewer plants were nodulated

330

when inoculated with bulk soil compared to WT (Fisher's exact test, P = 0.0076, P = 0.0686 and P <

331

0.0001 for P. sativum, P. vulgaris and G. soja, respectively; Figure 4E). Species identification of the

332

nodule isolates suggested that R. anhuiense is the microsymbiont of P. sativum, B. elkanii is the

333

microsymbiont of P. vulgaris, and B. japonicum is the microsymbiont of G. soja (Supplementary

334

Figure 13); among them R. anhuiense, and B. japonicum were AM symbiosis-promoted rhizobial

335

species in the M. truncatula rhizosphere (Figure 4B and Supplementary Figure 8). Collectively, our

336

results suggest that, in rhizobia-poor native soil, AM symbiosis-promoted rhizobia accumulation in

337

the rhizosphere might stimulate nodulation.

Jo

ur

na

lP

re

-p

ro

of

316

338

DISCUSSION Legume plants can form tripartite symbiotic associations with rhizobia and AM fungi. Here, we

340

studied the impact of AM and rhizobial symbiosis on the rhizosphere and root microbial

341

communities in M. truncatula grown in native soil using quantitative microbiota profiling (QMP).

342

The differentially expressed PSR genes and fatty acid biosynthetic genes between M+ and M- plants

343

are consistent with the role of AM fungi in facilitating the uptake of phosphorus, whereas fatty acids

344

synthesized in the host plants are transferred to the fungi to sustain mycorrhizal colonization (Jiang

345

et al., 2017; Luginbuehl et al., 2017). Genes related to the defense response and ROS (response to

346

oxidative stress) were more highly expressed in M- roots, suggesting an interaction between

347

mycorrhizas or phosphate stress and plant immunity (Castrillo et al., 2017). A weak, defective AM

348

phenotype has been observed in M. truncatula lyk3 and ipd3 mutants under controlled laboratory

349

conditions (Horvath et al., 2011; Jin et al., 2018; Zhang et al., 2015). Intriguingly, in native soil

350

conditions, our QMP analysis revealed that the rhizobial symbiotic mutants, lyk3, nfp and ipd3

351

exhibited increased AM colonization and phosphorus content in the roots compared to WT,

352

indicating that AM symbiosis might compensate for the loss of rhizobial symbiosis in M+R- plants in

353

native soil. Notably, this interesting phenomenon would not be detected by standard RMP analysis,

354

which does not consider the variability in microbial loads between samples. It will be important to

355

investigate competition and/or mutualist interactions between AM fungi and rhizobia in native soil in

356

the future.

Jo

ur

na

lP

re

-p

ro

of

339

357

Our quantitative analysis revealed that M+R- mutants have an elevated bacterial load in the

358

rhizophere accompanied by a reduced bacterial load in roots, whereas M-R- mutants have a reduced

359

bacterial load in the rhizosphere accompanied by an elevated bacterial load in roots. Niche

360

competition and nutrient exchange greatly influence fitness in multispecies interactions (Hibbing et

361

al., 2010). Our results suggest that AM fungi and the bacterial community may compete in the root

362

compartment, but nutrient transport via a hyphal network may favor bacterial colonization in the

363

rhizosphere compartment where the relatively larger space might alleviate competition. Moreover,

364

compared to WT, we observed that a broad-spectrum of bacterial phyla were consistently altered in

365

M+R- mutants or in M-R- mutants, indicating a robust association between the symbioses and

366

bacterial microbiota. Previous studies of the gut microbiota identified changes in microbial load as a

367

key driver of the altered microbiota observed in a cohort of patients with Crohn’s disease

368

(Vandeputte et al., 2017). Here we show that the altered bacterial load itself could be a key identifier

369

of a symbiosis-associated environment configuration. Based on RMP data, NFR5 (the L. japonicus ortholog of NFP) is required for normal root

371

microbiota assembly (Thiergart et al., 2019; Zgadzaj et al., 2016). Intriguingly, when we examined

372

the explained variance from a paired comparison of bacterial communities between WT and each of

373

the mutant lines determined by QMP, the M. truncatula NFP was not required to assemble a normal

374

root and rhizosphere microbiota. It is likely that microbiota composition will be different in soils

375

derived from different locations. Nevertheless, our results showed that the composition and diversity

376

of root-associated microbiota based on QMP and RMP matrices could be very different. Thus, it will

377

be critical to quantitatively track changes in the root-associated microbiota of wild type and mutant

378

genotypes of plants, and to link microbial load to quantitative features, such as physiological and

379

biochemical parameters, which cannot be addressed using relative abundance approaches (Guo et al.,

380

2019).

na

lP

re

-p

ro

of

370

In both leguminous and non-leguminous plants, species from the rhizobial genera were enriched

382

in the rhizosphere and root compartments relative to bulk soil samples (Hartman et al., 2017; Yeoh et

383

al., 2017). Comparative genomic analysis and studies with a tripartite plant system revealed that root

384

colonization for Rhizobiales isolates operate in a modular fashion and might be relevant to microbial

385

homeostasis in healthy roots (Garrido-Oter et al., 2018). Our statistical analyses revealed rhizobial

386

ASVs represent 45% of the consistently decreased ASVs in M- mutants and that the quantified

387

abundance of symbiotic rhizobial species was significantly lower in the rhizosphere soil of M- plants

388

compared to WT or M+ mutants, suggesting that AM symbiosis functions to structure rhizobial

389

communities. Moreover, Rhizobiales ASVs were more likely to be identified as microbial hubs in M+

390

plants compared to M- plants in both rhizosphere and root compartments. In this light, Rhizobiales

391

ASVs might mediate between the plant and the root-associated microbiota for the benefit of the plant

392

“holobiont” (Castrillo et al., 2017; Hassani et al., 2018; van der Heijden and Hartmann, 2016).

393

Jo

ur

381

AM fungi are important components of the soil ecosystem, extra-radical mycelium (ERM) extends

into the soil and acquires nutrients for the plant (van der Heijden et al., 2015). Our results show that

395

AM symbiosis promotes rhizobia accumulation in the rhizosphere, suggesting that the large surface

396

area of the ERM might provide a nutrient-rich niche to support the colonization and growth of

397

rhizobial communities. Our new findings are further supported by a previous study (Garrido-Oter et

398

al., 2018) which compared Rhizobiales among diverse plant hosts. Of these, AMF host plants Lotus

399

(12.43%) and barley (10.91%) displayed higher relative abundances of Rhizobiales in their

400

rhizospheres compared to the non-mycorrhizal plant Arabidopsis (4.45%), indicating that AM

401

symbiosis promotion of rhizobia accumulation in the rhizosphere might be a general mechanism to

402

further promote plant growth. The genetic program for AM has been recruited for legume-rhizobial

403

symbiosis (Oldroyd, 2013) and our data suggests that AM symbiosis also provides a

404

microenvironment, namely, the soil region influenced by the mycorrhizal roots which has been

405

termed as ‘mycorrhizosphere’ (Priyadharsini et al., 2016), for the co-evolution of AM fungi and

406

rhizobia.

lP

re

-p

ro

of

394

Jo

ur

na

407

408

METHODS

409

Field investigation of AM fungi in the botanical garden and soil sampling. At the start of this project, 27 samples from roots, rhizosphere soils and surrounding bulk soils

411

of 8 functionally distinct groups of naturally grown plants (Goosegrass, greenbean, mulberry, okra,

412

peanut, poplar, soja and soybean) were collected from a botanical garden in Nanjing Normal

413

University (N32°06'31.10", E118°54'21.81"). Through sequencing the 18S rRNA gene for AM fungi

414

of the phylum Glomeromycota, we identified 119 Glomeromycota 18S amplicon sequence variants

415

(ASVs) belonging to the 45 currently known AM fungal taxa [“virtual taxa” (VT)] in 21 of the

416

samples (data not shown). The physicochemical characteristics of the soils were the following: pH,

417

7.2; organic matter, 18.7 g/kg; available N, 60.6 mg/kg; available P, 6.1 mg/kg; available K, 116.7

418

mg/kg; soil texture, loam. Soils of this garden were sampled and transported back to the greenhouse

419

and mixed individually and homogenized thoroughly before use.

420

Plant material and plant growth conditions.

na

lP

re

-p

ro

of

410

The effects of AM symbiosis on microbial communities were surveyed using Medicago plants,

422

including wild type M. truncatula (‘Jemalong A17’), and the corresponding AM and

423

legume-rhizobial symbioses dysfunction mutants (lyk3-1, nfp-1, ipd3-1, dmi2-1, dmi3-1 and pt4-1).

424

LysM receptor kinase 3 (LYK3) and Nod factor perception (NFP) are Nod factor receptors required

425

for rhizobial symbiosis (Oldroyd, 2013). IPD3 is a member of the common symbiotic signaling

426

pathway required for rhizobial and mycorrhizal symbioses, however, the ipd3 mutant phenotype is

427

strongly dependent on the genetic background and ipd3-1 Jemalong plants displayed normal

428

arbuscule development (Oldroyd, 2013). Unlike lyk3 and nfp mutants, ipd3 was shown to be

429

impaired only in infection thread formation but not in root nodule organogenesis in M. truncatula,

430

and uninfected nodules (bump) were found on ipd3 roots (Jin et al., 2018). DMI2 is a symbiosis

431

receptor-like kinase, DMI3 is a calcium- and calmodulin dependent serine/threonine protein kinase,

432

both are required for rhizobial symbiosis and mycorrhizal establishment (Oldroyd, 2013). PT4 plays

433

a role in the acquisition of Pi released by the AM fungus and mutation of PT4 results in premature

Jo

ur

421

434

degradation of arbuscules in M. truncatula (Parniske, 2008). Seeds of the Medicago plants were mechanically scarified, surface-sterilized (10 s 75% EtOH, 2

436

min 25% NaClO and rinsed 5 times with sterile water) and germinated on sterile 1% water agar

437

plates in the dark (36 h at 4 °C and an additional 48 h at 28 °C). Then seedlings were transplanted

438

into pots containing soil batches described above. Two seedlings were transplanted into each pot and

439

12 pots were used for each plant genotype, 6 pots without plants were used as bulk soil controls. All

440

the pots with or without plants were randomly placed in greenhouse (16 h light /8 h dark period at

441

22ºC; photon flux density;250 μmol m2/s) and reset in the platform monthly. Plants were watered

442

with tap water according to their growth. The plants were harvested after 60 days, when all the plants

443

were still in their stem elongation developmental stage.

444

Collecting the rhizosphere and root samples of the Medicago plants.

re

-p

ro

of

435

To isolate rhizosphere and root microbial communities, for each pot, plant roots were cautiously

446

pulled out and separated from the adhering soil particles with sterilized forceps. The nodules were

447

removed prior to rhizosphere soil collection and root sonication. The roots were then replaced into 50

448

ml Falcon tubes containing 25 ml PBS buffer (140 mM NaCl, 2.7 mM KCl, 10 mM Na2HPO4, 1.8

449

mM KH2PO4) to obtain rhizosphere and root samples as previously described (Edwards et al., 2015;

450

Wang et al., 2020). That is, roots in the first wash step were stirred vigorously to clean all of the soil

451

from the root surfaces, this soil was then collected after centrifugation (20 min at 4,000 rpm) and

452

treated as rhizosphere soil. Root samples were obtained after 3 sonication procedures, for each

453

procedure, the cleaned roots were placed in fresh PBS buffer for sonication (30 s at 42 Hz, 30 s pause)

454

to remove tightly adhering soil particles. For each plant genotype, 5 biological repeats, which were

455

defined as 5 individual pots randomly picked among the 12 replicates, were processed.

Jo

ur

na

lP

445

456

To samples collected for gene expression analysis (RNA-seq), plant roots were cautiously

457

pulled out and washed in PBS buffer, then nodules were removed and the remaining clean roots were

458

promptly frozen in liquid nitrogen. Frozen roots were stored at -80℃ until RNA extraction. For each

459

plant genotype, 5 biological repeats, which was defined as 5 individual pots randomly picked among

460

the 12 replicates, were processed for RNA extraction.

461

Design of bacteria 16S, rhizobia rpoB and Glomeromycota 18S synthetic chimeric DNA spikes. The methods for designing synthetic chimeric DNA spikes were previously described (Tkacz et

463

al., 2018; Wang et al., 2020). That is, bacteria 16S synthetic spike was designed based on the

464

prokaryotic 16S rRNA gene V3–V4 regions primers 341F-806R (Roggenbuck et al., 2014), with a

465

stuffer sequence of 429 bp matched to the length of its natural PCR products. The single copy rpoB

466

gene primers rpoB1479F (5ʹ -GAT CGA AAC GCC GGA AGG-3ʹ) and rpoB1831R (5ʹ-TGC ATG

467

TTC GAG CCC AT-3ʹ), which used for the intraspecific level rhizobial diversity analysis in our

468

previous studies (Wang et al., 2018; Zhang et al., 2017), was used to design rhizobia rpoB synthetic

469

spike, with a stuffer sequence of 343 bp matched to the length of its natural PCR products.

470

Glomeromycota 18S rRNA gene primers NS31 (5ʹ-TTG GAG GGC AAG TCT GGT GCC-3ʹ) and

471

AML2 (5ʹ-GAA CCC AAA CAC TTT GGT TTC C-3ʹ), which was used for assessment of AM

472

fungal diversity (Davison et al., 2015), was used to design the Glomeromycota 18S synthetic spike,

473

with a stuffer sequence of 524 bp matched to the length of its natural PCR products. All of the

474

synthetic spikes were concordantly integrated into a synthetic chimeric DNA sequence (1541 bp) and

475

synthesised by Sangong Biotech. Plasmid pUC57 carrying our synthetic spikes was transformed into

476

Escherichia coli strain TOP10. The recombinant plasmids in TOP10 were extracted, purified and

477

with number of copies quantified according to its concentration and molecular weight (2.15 × 105

478

pg-1), then diluted to an appropriate concentration before use.

479

Sample quantitation, DNA extraction and sequencing.

Jo

ur

na

lP

re

-p

ro

of

462

480

A sample quantitation method was previously described (Wang et al., 2020). Namely, bulk soil,

481

rhizosphere and root samples were weighed, and the synthetic chimeric DNA spikes were added into

482

the samples before DNA extraction. For whole bacteria and rhizobia surveys, 500 pg (about 108

483

million spike copies) of synthetic spikes were added to both root and rhizosphere samples prior to

484

DNA isolation. For the Glomeromycota survey, 5 pg of synthetic spike was added to root samples

485

prior to DNA isolation. The amounts of spikes added for bacteria and fungi quantitation were

486

different based on the microbial loads of our samples (Wang et al., 2020) and suggestions from

487

Tkacz et al (Tkacz et al., 2018). DNA extraction was performed using FastDNA SPIN Kit for Soil

488

(MP Biomedicals). The V3–V4 regions of bacterial 16S rRNA gene sequences were amplified from

489

DNA extracts using the 341F-806R universal primers (Roggenbuck et al., 2014), rhizobial rpoB gene

490

sequences were amplified using the rpoB1479F-rpoB1831R universal primers (Zhang et al., 2017),

491

and Glomeromycota nuclear 18S rRNA gene sequences were amplified using the NS31-AML2

492

universal primers (Davison et al., 2015). The library construction and sequencing (Illumina MiSeq,

493

PE 2 × 300 bp) was conducted by the commercial service provider Shanghai Hanyu Biotech lab

494

(Shanghai, China). For the bacterial 16S rRNA gene, a total of 2,336,892 high-quality sequences were obtained

496

with a median read count per sample of 20,810 (range: 6,017–56,633). For the bacterial rpoB gene, a

497

total of 2,472,269 high-quality sequences were obtained with a median read count per sample of

498

33,421 (range: 13,897–90,544). For Glomeromycota 18S rRNA gene, a total of 1,837,696

499

high-quality sequences were obtained with a median read count per sample of 49,041 (range:

500

20,920–102,529). Raw microbiota sequencing fastq data reported in this study are publicly deposited

501

in the SRA (Sequence Read Archive) database under the BioProject PRJNA551144 (bacteria 16S

502

rRNA gene accession no. SAMN12141687-SAMN12141761, rhizobia rpoB gene accession no.

503

SAMN12141762-SAMN12141836

504

SAMN12141837-SAMN12141871).

505

RNA extraction and sequencing.

na

lP

re

-p

ro

of

495

Glomeromycota

18S

rRNA

gene

accession

no.

Jo

ur

and

506

To extract RNA for gene expression analysis, the complete frozen root systems of plants were

507

pulverized in liquid nitrogen and RNAs were extracted using Trizol® Reagent. The quality and

508

integrity of RNAs were determined by agarose gel electrophoresis and a 2100 Bioanalyzer. The

509

qualified RNAs were treated with DNase (5 U/μL) (TaKaRa, Japan) at 37℃ for 30 mins. The DNase

510

treated RNAs were purified by Dynabeads® Oligo (dT)25 (Life，USA). The cDNA library

511

construction and sequencing (Illumina HiSeq, PE 2 × 150 bp) was performed by the commercial

512

service provider Shanghai Hanyu Biotech lab at Shanghai (Shanghai, China). For each of 5 replicate

513

plants per genotype, we generated 23,977,894 ± 6,809,822 clean reads (total = 839,226,320 reads).

514

Raw RNA sequencing fastq data reported in this study are available under the BioProject

515

PRJNA551144, accession no. SAMN12233320-SAMN1223354.

516

Processing of amplicon data. Raw Illumina fastq files were quality filtered and taxonomy analyzed using QIIME 2 Core

518

2018.11 distribution (Bolyen et al., 2019). For bacterial 16S rRNA gene amplicon sequences, all

519

unassigned sequences and sequences annotated as mitochondria, chloroplast, Archaea and

520

unclassified were removed in downstream analysis. For rpoB amplicon sequences, only sequences

521

belonging to the clade of Proteobacteria were retained in downstream analysis, then rhizobia (both

522

symbiotic and non-symbiotic rhizobia) reads at the genus level and potential symbiotic rhizobia

523

reads at the intraspecific level were counted and annotated. The species-level positions of rpoB ASVs

524

were annotated according to the rpoB references (Wang et al., 2020) by using the

525

classify-consensus-blast method (--p-maxaccepts 1000, --p-perc-identity 0.977, --p-query-cov 0.8

526

and --p-min-consensus 0.8), as 97.7% of rpoB sequence similarity correlated with membership of

527

different species (with a DDH value < 70% and an ANI value < 94.3 %) (Adekambi et al., 2008; Vos

528

et al., 2012). For Glomeromycota 18S rRNA processing, Glomeromycota NS31-AML2 amplicon

529

sequences were BLASTed against our Glomeromycota 18S spike sequence and against

530

Glomeromycota 18S rRNA gene sequences in the MaarjAM database (Opik et al., 2010) using an

531

approach similar to those described by (Davison et al., 2015). Only sequences with a best hit with ≥

532

97% similarity against an AM virtual taxon (VT) in the MaarjAM database and an alignment length

533

≥ 280 bp were counted and annotated. After this processing, the filtered representative sequences and

534

ASV tables were then used to determine taxonomic abundances and subsequent statistical analyses in

535

R.

536

Quantitation and copy number weighting of ASV tables.

Jo

ur

na

lP

re

-p

ro

of

517

537

To quantify microbial abundance in rhizosphere and root samples, firstly, the quantified

538

abundance of genes (16S rRNA, rpoB and 18S rRNA) for each ASV per gram of sample (QAGA)

539

was estimated using the following equation: QAGA =

S × SA SS

540

where i represents the IDs of each ASV, S is the number of synthetic spikes added per gram of

541

sample. SS is the number of sequenced reads assigned to synthetic spikes, SA is the number of

542

sequenced reads assigned to each ASV. For the multi-copy 16S rRNA gene, the 16S rRNA gene copy number of each ASV (CNAi) was

544

estimated using rrnDB database (Stoddard et al., 2015). In brief, the mean gene copy number (if

545

available) of the immediate child taxa with a bootstrap confidence estimate above the 80% threshold

546

for each ASV was used as the mean copy number for that ASV. For any ASV without copy number

547

data, the mean copy number of its parent was used for that ASV. At present, the RDP classifier also

548

provides gene copy number adjustment for 16S rRNA gene sequences to better estimate bacterial

549

relative abundance (Wang et al., 2007). Finally, the quantified abundance of bacteria for each ASV

550

(QABA) and the quantified abundance of total bacteria (QAB) per gram of sample was estimated

551

using the following equations: QAGA , CNA

lP

QABA =

re

-p

ro

of

543

QAB =

QABA

For the single copy rpoB gene, QAGAs per gram of sample are equal to QABAs per gram of

553

sample. The quantified abundance of rhizobia (genus of Rhizobium, Bradyrhizobium, Mesorhizobium

554

and Sinorhizobium/Ensifer) in our samples were compared between single copy rpoB gene and 16S

555

rRNA gene surveys to evaluate quantitative accuracy (Wang et al., 2020). Glomeromycota 18S ASVs

556

were quantified but the copy number was not weighted due to the lack of a corresponding database.

557

Processing of RNA-seq data.

Jo

ur

na

552

558

Illumina raw sequencing reads from RNA samples were preprocessed using Trimmomatic v0.35

559

(Bolger et al., 2014), reads filtered with a quality cutoff of 20 and shorter than 75 bp were discarded.

560

Transcript-level expression analysis of the clean RNA reads was performed using HISAT and

561

StringTie (Pertea et al., 2016). Namely, clean reads were mapped to the genome of M. truncatula

562

A17 (Mt4.0; www.medicagohapmap.org) by HISAT2 software. The expression level of each gene

563

mapped to A17 was calculated by StringTie software, then Cuffdiff software was used to calculate

564

differentially expressed genes (DEGs) between plant genotypes. Genes with value of |log2fold

565

change

| > 1.5 and q-value < 0.01 was defined as the DEGs between plant genotypes in this study.

566

Statistical Analysis For RNA-seq data, the hierarchical clustering analyses were performed with the ‘Heatmap’

568

function in R from the ComplexHeatmap package (Gu et al., 2016), using the sets of DEGs identified

569

in our experiment. Genes were clustered on the basis of the Euclidean distance and with the

570

complete-linkage method. Genes belonging to each cluster were submitted to Gene Ontology (GO)

571

enrichment analyses on the PlantGSEA platform (Yi et al., 2013) to identify over-represented

572

biological processes. For microbiota data, CAP (canonical analysis of principal coordinates)

573

(Anderson and Willis, 2003) plots based on the Bray-Curtis dissimilarities of quantified bacteria

574

were used to visualize genotype relationships. ANOSIM (analysis of similarities) (CLARKE, 1993)

575

with 9999 permutations were used to test significant differences between groups based on the Bray–

576

Curtis dissimilarities. Differential taxa with relative/absolate abundance between wild type A17 and

577

symbiosis mutants were identified via the nonparametric Kruskal–Wallis rank test with false

578

discovery rate correction (α= 0.05), using the “kruskal” function in the package “agricolae” in R

579

software (Mendiburu, 2010). The significance taxa with false discovery rate (FDR) <0.05 were

580

assigned

581

enriched/depleted taxa in symbiosis mutants. Microbial co-abundance network analyses were

582

performed using the igraph R package (Csardi and Nepusz, 2005), namely, the significant correlation

583

(Spearman's correlation test, PFDR < 0.05) between paired ASVs based on the QMP data set of 16S

584

rRNA genes were included into the network analyses, and the topological parameters of the network,

585

including the degree, betweeness centrality and closeness centrality were calculated. Microbial hubs

586

were identified according to the method provided by Agler et. al. (Agler et al., 2016). In brief, a

587

normal distribution of the log-transformed degree, betweenness centrality or closeness centrality data

588

was used to calculate the values above which nodes (ASVs) can be considered outliers corresponding

589

to p < 0.2. Nodes (ASVs) that were above this value for all three parameters were considered to be

590

hubs. Network visualizations were constructed using Gephi (Bastian et al., 2009). Data visualization

591

was primarily generated using ggplot2 R package (Wickham, 2016). Heatmaps were generated using

592

the ComplexHeatmap R package (Gu et al., 2016). Annotated phylogenetic trees were carried out

593

using online EvolView (Zhang et al., 2012).

either

absolute

abundance

elevated/decreased

taxa

or

relative

abundance

Jo

ur

as

na

lP

re

-p

ro

of

567

594

Rhizosphere inoculation experiments and species identification of the nodule isolates Four legumes, namely Vicia faba, Pisum sativum, Phaseolus vulgaris and Glycine soja, were

596

used for the rhizosphere inoculation experiment. Seed of the legumes were surface sterilized and

597

germinated as described above for M. truncatula, and inoculated with 1 ml 105 diluted bulk soil and

598

rhizosphere microbial communities from M. truncatula WT and mutant plants (lyk3, dmi3 and pt4).

599

Plants were grown in sterilized sand (sand/perlite = 3:1, sterilized by autoclave on two consecutive

600

days for 20 min at 121°C) in Leonard jars (Vincent, 1970) in greenhouse. At 21 dpi, plants were

601

harvested and the nodule numbers were determined. To identify the root nodule microsymbionts of

602

these legumes, root nodules were crushed, spread onto TY solid medium (tryptone at 5 g/L, yeast

603

extract at 3 g/L, and CaCl2 at 0.6 g/L) and incubated at 28 °C for 5 days. The nodule isolates were

604

then purified twice and identified by rpoB sequencing, as described earlier (Wang et al., 2018).

605

Data Availability

ro

-p

re lP

Raw microbiota and RNA sequencing fastq data reported in this study are publicly deposited in

na

606

of

595

the SRA database under the BioProject PRJNA551144 as described above. The R scripts, rpoB

608

reference database, as well as the metadata of bacterial 16S rRNA gene, rhizobial rpoB gene and

609

Glomeromycota 18S rRNA gene employed in the microbiota analyses are available at

610

https://github.com/godlovexiaolin/Mycorrhizal-symbiosis-modulates-the-rhizosphere-microbiota-to-

611

promote-rhizobia-legume-symbiosis.

Jo

ur

607

612 613

FUNDING

614

The research was supported by Chinese Academy of Sciences (ZDRW-ZS-2019-2), National

615

Science Foundation (31825003, 31730103 and 31970323), The Strategic Priority Research Program

616

“Molecular Mechanism of Plant Growth and Development” of the Chinese Academy of Sciences

617

(XDB27040207) and China National GeneBank (CNGB).

618

619

AUTHOR CONTRIBUTIONS

620

Conceptualization, X. W. and E. W.; Methodology, X. W.; Investigation, X. W., H. F., L. W., Y.

621

W., M. W., X. X., H. C., J. Q. and K. S.; Writing – Original Draft, X. W. and E. W.; Writing –Review

622

& Editing, X. W. and E. W.; Funding Acquisition, H. L., X. Z. and E. W.; Resources, W. H., C. W., C.

623

D., Z. C., C. T. and N. Y.; Supervision, H. L., X. Z. and E. W.

624

of

626

DECLARATION OF INTERESTS The authors declare no competing interests.

ro

625

Jo

ur

na

lP

re

-p

627

628

REFERENCE

629

Adekambi, T., Shinnick, T.M., Raoult, D., and Drancourt, M. (2008). Complete rpoB gene

630

sequencing as a suitable supplement to DNA-DNA hybridization for bacterial species and genus

631

delineation. Int. J. Syst. Evol. Microbiol. 58:1807-1814. Afkhami, M.E., and Stinchcombe, J.R. (2016). Multiple mutualist effects on genomewide

633

expression in the tripartite association between Medicago truncatula, nitrogen-fixing bacteria

634

and mycorrhizal fungi. Mol. Ecol. 25:4946-4962.

of

632

Agler, M.T., Ruhe, J., Kroll, S., Morhenn, C., Kim, S.T., Weigel, D., and Kemen, E.M. (2016).

636

Microbial hub taxa link host and abiotic factors to plant microbiome variation. PLoS. Biol.

637

14:e1002352.

-p

re

639

Anderson, M.J., and Willis, T.J. (2003). Canonical analysis of principal coordinates: A useful method of constrained ordination for ecology. Ecology 84:511-525.

lP

638

ro

635

Ballesteros-Almanza, L., Altamirano-Hernandez, J., Peña Cabriales, J., Santoyo, G.,

641

Sanchez-Yañez, J., Valencia-Cantero, E., Macias-Rodriguez, L., López-Bucio, J., Cárdenas

642

Navarro, R., and Farías-Rodríguez, R. (2010). Effect of co-inoculation with mycorrhiza and

643

rhizobia on the nodule trehalose content of different bean genotypes. MicrobiologyOpen

644

4:83-92.

646 647 648

ur

Jo

645

na

640

Bastian, M., Heymann, S., and Jacomy, M. (2009). Gephi: An open source software for exploring and manipulating networks. Bolger, A.M., Lohse, M., and Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30:2114-2120.

649

Bolyen, E., and Rideout, J.R., and Dillon, M.R., and Bokulich, N.A., and Abnet, C.C., and

650

Al-Ghalith, G.A., and Alexander, H., and Alm, E.J., and Arumugam, M., and Asnicar, F.,

651

et al. (2019). Reproducible, interactive, scalable and extensible microbiome data science using

652

QIIME 2. Nat. Biotechnol. 37:852-857.

653 654

Bonfante, P., and Anca, I.A. (2009). Plants, mycorrhizal fungi, and bacteria: A network of interactions. Annu. Rev. Microbiol. 63:363-383.

655

Bulgarelli, D., Schlaeppi, K., Spaepen, S., Ver Loren van Themaat, E., and Schulze-Lefert, P.

656

(2013). Structure and functions of the bacterial microbiota of plants. Annu. Rev. Plant Biol.

657

64:807-838.

658

Carrion, V.J., Perez-Jaramillo, J., Cordovez, V., Tracanna, V., de Hollander, M., Ruiz-Buck, D.,

659

Mendes, L.W., van Ijcken, W.F.J., Gomez-Exposito, R., Elsayed, S.S., et al. (2019).

660

Pathogen-induced activation of disease-suppressive functions in the endophytic root

661

microbiome. Science 366:606-612. Case, R.J., Boucher, Y., Dahllof, I., Holmstrom, C., Doolittle, W.F., and Kjelleberg, S. (2007).

663

Use of 16S rRNA and rpoB genes as molecular markers for microbial ecology studies. Appl.

664

Environ. Microbiol. 73:278-288.

ro

of

662

Castrillo, G., Teixeira, P.J., Paredes, S.H., Law, T.F., de Lorenzo, L., Feltcher, M.E., Finkel,

666

O.M., Breakfield, N.W., Mieczkowski, P., Jones, C.D., et al. (2017). Root microbiota drive

667

direct integration of phosphate stress and immunity. Nature 543:513-518.

lP

re

-p

665

Catford, J.G., Staehelin, C., Lerat, S., Piche, Y., and Vierheilig, H. (2003). Suppression of

669

arbuscular mycorrhizal colonization and nodulation in split-root systems of alfalfa after

670

pre-inoculation and treatment with Nod factors. J. Exp. Bot. 54:1481-1487.

672 673 674

ur

Clarke, K.R. (1993). Non-parametric multivariate analyses of changes in community structure. Aust.

Jo

671

na

668

Journal Ecol. 18:117-143.

Csardi, G., and Nepusz, T. (2005). The igraph software package for complex network research. InterJournal Complex Systems 1695.

675

Davison, J., Moora, M., Opik, M., Adholeya, A., Ainsaar, L., Ba, A., Burla, S., Diedhiou, A.G.,

676

Hiiesalu, I., Jairus, T., et al. (2015). Global assessment of arbuscular mycorrhizal fungus

677

diversity reveals very low endemism. Science 349:970-973.

678

Edwards, J., Johnson, C., Santos-Medellin, C., Lurie, E., Podishetty, N.K., Bhatnagar, S., Eisen,

679

J.A., and Sundaresan, V. (2015). Structure, variation, and assembly of the root-associated

680

microbiomes of rice. Proc. Natl. Acad. Sci. U S A 112:E911-920.

681 682

Faust, K., and Raes, J. (2012). Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10:538-550.

683

Garrido-Oter, R., Nakano, R.T., Dombrowski, N., Ma, K.W., AgBiome, T., McHardy, A.C., and

684

Schulze-Lefert, P. (2018). Modular traits of the Rhizobiales root microbiota and their

685

evolutionary relationship with symbiotic rhizobia. Cell Host Microbe 24:155-167.

686

Gobbato, E., Marsh, John F., Vernié, T., Wang, E., Maillet, F., Kim, J., Miller, J.B., Sun, J.,

687

Bano, S.A., Ratet, P., et al. (2012). A GRAS-type transcription factor with a specific function

688

in mycorrhizal signaling. Curr. Biol. 22:2236-2241.

690

Gu, Z.G., Eils, R., and Schlesner, M. (2016). Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32:2847-2849.

of

689

Guo, X., Zhang, X., Qin, Y., Liu, Y.-X., Zhang, J., Zhang, N., Wu, K., Qu, B., He, Z., Wang, X.,

692

et al. (2019). Host-associated quantitative abundance profiling reveals the microbial load

693

variation of root microbiome. Plant Commun. 1:100003.

-p

ro

691

Hao, Z., Xie, W., Jiang, X., Wu, Z., Zhang, X., and Chen, B. (2019). Arbuscular mycorrhizal

695

fungus improves rhizobium–glycyrrhiza seedling symbiosis under drought stress. Agronomy

696

9:572.

lP

re

694

Hartman, K., van der Heijden, M.G., Roussely-Provent, V., Walser, J.C., and Schlaeppi, K.

698

(2017). Deciphering composition and function of the root microbiome of a legume plant.

699

Microbiome 5:2.

701 702 703

ur

Jo

700

na

697

Hassani, M.A., Duran, P., and Hacquard, S. (2018). Microbial interactions within the plant holobiont. Microbiome 6:58. Hibbing, M.E., Fuqua, C., Parsek, M.R., and Peterson, S.B. (2010). Bacterial competition: Surviving and thriving in the microbial jungle. Nat. Rev. Microbiol. 8:15-25.

704

Horvath, B., Yeun, L.H., Domonkos, A., Halasz, G., Gobbato, E., Ayaydin, F., Miro, K., Hirsch,

705

S., Sun, J.H., Tadege, M., et al. (2011). Medicago truncatula IPD3 is a member of the common

706

symbiotic signaling pathway required for rhizobial and mycorrhizal symbioses. Mol.

707

Plant-Microbe Interact. 24:1345-1358.

708

Javot, H., Penmetsa, R.V., Terzaghi, N., Cook, D.R., and Harrison, M.J. (2007). A Medicago

709

truncatula phosphate transporter indispensable for the arbuscular mycorrhizal symbiosis. Proc.

710

Natl. Acad. Sci. U S A 104:1720-1725.

711

Jiang, Y., Wang, W., Xie, Q., Liu, N., Liu, L., Wang, D., Zhang, X., Yang, C., Chen, X., Tang, D.,

712

et al. (2017). Plants transfer lipids to sustain colonization by mutualistic mycorrhizal and

713

parasitic fungi. Science 356:1172-1175.

714

Jin, Y., Chen, Z.X., Yang, J., Mysore, K.S., Wen, J.Q., Huang, J.R., Yu, N., and Wang, E.T.

715

(2018). IPD3 and IPD3L function redundantly in rhizobial and mycorrhizal symbioses. Front.

716

Plant Sci. 9:267. Li, F.N., Liu, S.W., Lu, Q.P., Zheng, H.Y., Osterman, I.A., Lukyanov, D.A., Sergiev, P.V.,

718

Dontsova, O.A., Liu, S.S., Ye, J.J., et al. (2019a). Studies on antibacterial activity and

719

diversity of cultivable actinobacteria isolated from mangrove soil in Futian and Maoweihai of

720

China. Evid-Based Compl. Alt.

ro

of

717

Li, N., Lin, B., Wang, H., Li, X., Yang, F., Ding, X., Yan, J., and Chu, Z. (2019b). Natural

722

variation in ZmFBL41 confers banded leaf and sheath blight resistance in maize. Nature Genet.

723

51:1540-1548.

lP

re

-p

721

Lindsay, P.L., Williams, B.N., MacLean, A., and Harrison, M.J. (2019). A phosphate-dependent

725

requirement for transcription factors IPD3 and IPD3L during arbuscular mycorrhizal symbiosis

726

in Medicago truncatula. Mol. Plant-Microbe Interact. 32:1277-1290.

ur

na

724

Luginbuehl, L.H., Menard, G.N., Kurup, S., Van Erp, H., Radhakrishnan, G.V., Breakspear, A.,

728

Oldroyd, G.E.D., and Eastmond, P.J. (2017). Fatty acids in arbuscular mycorrhizal fungi are

729

synthesized by the host plant. Science 356:1175-1178.

730 731 732 733 734 735 736 737 738

Jo

727

Martin, F.M., Uroz, S., and Barker, D.G. (2017). Ancestral alliances: Plant mutualistic symbioses with fungi and bacteria. Science 356:eaad4501. Mendiburu, F. (2010). Agricolae: Statistical procedures for agricultural research. R package version 1:1-8. Miransari, M. (2011). Interactions between arbuscular mycorrhizal fungi and soil bacteria. Appl. Microbiol. Biotechnol. 89:917-930. Oldroyd, G.E. (2013). Speak, friend, and enter: Signalling systems that promote beneficial symbiotic associations in plants. Nat. Rev. Microbiol. 11:252-263. Opik, M., Vanatoa, A., Vanatoa, E., Moora, M., Davison, J., Kalwij, J.M., Reier, U., and Zobel,

739

M. (2010). The online database MaarjAM reveals global and ecosystemic distribution patterns

740

in arbuscular mycorrhizal fungi (Glomeromycota). New Phytol. 188:223-241.

741 742

Parniske, M. (2008). Arbuscular mycorrhiza: The mother of plant root endosymbioses. Nat. Rev. Microbiol. 6:763-775.

743

Pertea, M., Kim, D., Pertea, G.M., Leek, J.T., and Salzberg, S.L. (2016). Transcript-level

744

expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc.

745

11:1650. Priyadharsini, P., Rojamala, K., Ravi, R.K., Muthuraja, R., Nagaraj, K., and Muthukumar, T.

747

(2016). Mycorrhizosphere: The extended rhizosphere and its significance. In: Plant-Microbe

748

Interaction: An approach to sustainable agriculture--Choudhary, D.K., Varma, A., and Tuteja, N.,

749

eds. Singapore: Springer Singapore. 97-124.

-p

ro

of

746

Roggenbuck, M., Baerholm Schnell, I., Blom, N., Baelum, J., Bertelsen, M.F., Sicheritz-Ponten,

751

T., Sorensen, S.J., Gilbert, M.T., Graves, G.R., and Hansen, L.H. (2014). The microbiome

752

of New World vultures. Nat. Commun. 5:5498.

lP

re

750

Shrihari, P., Sakamoto, K., Inubushi, K., and Akao, S. (2000). Interaction between

754

supernodulating or non-nodulating mutants of soybean and two arbuscular mycorrhizal fungi.

755

Mycorrhiza 10:101-106.

Jo

ur

na

753

756

Stoddard, S.F., Smith, B.J., Hein, R., Roller, B.R.K., and Schmidt, T.M. (2015). rrnDB:

757

Improved tools for interpreting rRNA gene abundance in bacteria and archaea and a new

758

foundation for future development. Nucleic Acids Res. 43:D593-D598.

759

Thiergart, T., Zgadzaj, R., Bozsoki, Z., Garrido-Oter, R., Radutoiu, S., and Schulze-Lefert, P.

760

(2019). Lotus japonicus symbiosis genes impact microbial interactions between symbionts and

761

multikingdom commensal communities. mBio 10: e01833-19.

762 763

Tkacz, A., Hortala, M., and Poole, P.S. (2018). Absolute quantitation of microbiota abundance in environmental samples. Microbiome 6:110.

764

van der Heijden, M.G., de Bruin, S., Luckerhoff, L., van Logtestijn, R.S., and Schlaeppi, K.

765

(2016). A widespread plant-fungal-bacterial symbiosis promotes plant biodiversity, plant

766

nutrition and seedling recruitment. ISME J 10:389-399.

767 768 769 770

van der Heijden, M.G., and Hartmann, M. (2016). Networking in the plant microbiome. PLoS. Biol. 14:e1002378. van der Heijden, M.G.A., Martin, F.M., Selosse, M.A., and Sanders, I.R. (2015). Mycorrhizal ecology and evolution: The past, the present, and the future. New Phytol. 205:1406-1423.

771

Vandeputte, D., Kathagen, G., D'Hoe, K., Vieira-Silva, S., Valles-Colomer, M., Sabino, J., Wang,

772

J., Tito, R.Y., De Commer, L., Darzi, Y., et al. (2017). Quantitative microbiome profiling links

773

gut community variation to microbial load. Nature 551:507-511.

775

Vincent, J.M. (1970). A manual for the practical study of the root-nodule bacteria: IBP Handbk 15

of

774

Oxford and Edinburgh: Blackwell Scientific Publications.

Vos, M., Quince, C., Pijl, A.S., de Hollander, M., and Kowalchuk, G.A. (2012). A comparison of

777

rpoB and 16S rRNA as markers in pyrosequencing studies of bacterial diversity. PloS One

778

7:e30600.

re

-p

ro

776

Wang, Q., Garrity, G.M., Tiedje, J.M., and Cole, J.R. (2007). Naive Bayesian classifier for rapid

780

assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol.

781

73:5261-5267.

na

lP

779

Wang, X., Wang, M., Xie, X., Guo, S., Zhou, Y., Zhang, X., Yu, N., and Wang, E. (2020). An

783

amplification-selection model for quantified rhizosphere microbiota assembly. Sci. Bull.

784

65:983-986

Jo

ur

782

785

Wang, X.L., Cui, W.J., Feng, X.Y., Zhong, Z.M., Li, Y., Chen, W.X., Chen, W.F., Shao, X.M.,

786

and Tian, C.F. (2018). Rhizobia inhabiting nodules and rhizosphere soils of alfalfa: A strong

787

selection of facultative microsymbionts. Soil Biol. Biochem. 116:340-350.

788

Wang, X.R., Pan, Q.A., Chen, F.X., Yan, X.L., and Liao, H. (2011). Effects of co-inoculation with

789

arbuscular mycorrhizal fungi and rhizobia on soybean growth as related to root architecture and

790

availability of N and P. Mycorrhiza 21:173-181.

791

Wickham, H. (2016). ggplot2: Elegant graphics for data analysis: Springer-Verlag New York.

792

Xue, L., Almasio, J., Fabianska, I., Saridis, G., and Bucher, M. (2019). Dysfunction in the

793

arbuscular mycorrhizal symbiosis has consistent but small effects on the establishment of the

794

fungal microbiota in Lotus japonicus. New Phytol. 224:409-420.

795

Yeoh, Y.K., Dennis, P.G., Paungfoo-Lonhienne, C., Weber, L., Brackin, R., Ragan, M.A.,

796

Schmidt, S., and Hugenholtz, P. (2017). Evolutionary conservation of a core root microbiome

797

across plant phyla along a tropical soil chronosequence. Nat. Commun. 8.

798 799

Yi, X., Du, Z., and Su, Z. (2013). PlantGSEA: A gene set enrichment analysis toolkit for plant community. Nucleic Acids Res. 41:W98-W103. Zgadzaj, R., Garrido-Oter, R., Jensen, D.B., Koprivova, A., Schulze-Lefert, P., and Radutoiu, S.

801

(2016). Root nodule symbiosis in Lotus japonicus drives the establishment of distinctive

802

rhizosphere, root, and nodule bacterial communities. Proc. Natl. Acad. Sci. U S A

803

113:E7996-E8005.

of

800

Zhang, H., Gao, S., Lercher, M.J., Hu, S., and Chen, W.H. (2012). EvolView, an online tool for

805

visualizing, annotating and managing phylogenetic trees. Nucleic Acids Res. 40:W569-572.

806

Zhang, X.W., Dong, W.T., Sun, J.H., Feng, F., Deng, Y.W., He, Z.H., Oldroyd, G.E.D., and

807

Wang, E.T. (2015). The receptor kinase CERK1 has dual functions in symbiosis and immunity

808

signalling. Plant J. 81:258-267.

lP

re

-p

ro

804

Zhang, X.X., Guo, H.J., Jiao, J., Zhang, P., Xiong, H.Y., Chen, W.X., and Tian, C.F. (2017).

810

Pyrosequencing of rpoB uncovers a significant biogeographical pattern of rhizobial species in

811

soybean rhizosphere. J Biogeogr. 44:1491-1499.

ur

Jo

812

na

809

FIGURES LEGENDS

814

Figure 1. Microbiota and Transcriptome Profiles of Wild Type A17, Rhizobial and Arbuscular

815

Mycorrhizal Symbiosis Defective Mutants Grown in Native soils.

816

(A and B) Estimated the microbial load of Glomeromycota 18S rRNA genes (A) and the total

817

phosphorus content (B) in roots of different plant genotypes. Statistical significance was determined

818

via multiple comparison with Kruskal-Wallis test and with false discovery rate correction (α = 0.05).

819

Post hoc test is indicated by letters at the top, sample groups with the same letter are

820

indistinguishable at 95% confidence. Data are mean ± s.e.m., n = 5 biological replicates.

821

(C) Hierarchical clustering of 3,606 genes (DEGs) that were differentially expressed in roots of M.

822

truncatula genotypes. The optimal number of clusters (7 clusters) was determined via k-means

823

partitioning. Columns on the right indicate genes that are categorized as symbiosis affected gene sets

824

(rhizobia only affected genes, mycorrhizal fungi only affected genes, additive and non-additive

825

affected genes) in non-natural condition using the tripartite mutualism system (Afkhami and

826

Stinchcombe, 2016).

827

(D and E) Relative versus quantitative microbiota profiling. (D) Relative bacterial microbiota

828

profiles deduced from standard microbiota sequencing protocols. (E) Quantitative bacterial

829

microbiota profiles deduced from complementing sequencing with microbial cell counts (cells per

830

gram of samples). The top 8 most abundant phyla are depicted, with all others pooled into ‘Other’.

831

Proteobacteria are shown at the class level due to its high load. The total bacterial loads are shown

832

above the bar charts for each sample category. Pound sign (#) indicates bacterial phyla that are

833

significantly higher in mutants compared to WT at P < 0.05; asterisk (*) indicates bacterial phyla that

834

are significantly lower in mutants vs WT at P < 0.05.

835

Figure 2. AM Symbiosis is Required to Assemble a Normal Root-Associated Microbiota in

836

Native Soil.

837

(A and B) CAP biplots showing genotype variation of quantified bacterial taxa in rhizosphere (A)

838

and root samples (B) using Bray-Curtis distances constrained by plant genotype. Plant genotypes that

Jo

ur

na

lP

re

-p

ro

of

813

can form proficient or defective AM symbiosis are indicated with triangles and circles, respectively.

840

The confidence ellipses include 68% of samples for each genotype.

841

(C) Heatmap and phylogenetic tree of the 85 bacterial ASVs with microbial load increased or

842

decreased in rhizosphere soils of symbiosis mutants compared to those of wild type (WT). The ML

843

phylogenetic tree of bacterial ASVs was constructed using ASV representative sequences (1000

844

bootstrap replications), and tree branches were colored according to its taxonomy at the phylum level.

845

Color key of heatmap shown as (–)log10FDR represents the increase (red) and decrease (dark green)

846

of ASV load in rhizosphere soils of mutants compared to WT. The asterisk represent ASVs

847

consistently increased or decreased in rhizosphere soils of M- mutants compared with wild type

848

(WT).

849

Figure 3. AM Symbiosis Impacts Rhizobiales-Hubs in the Quantitative Microbial

850

Co-Abundance Networks.

851

(A and B) Co-abundance network of Spearman’s rho correlation between paired ASVs based on the

852

QMP data set of 16S rRNA genes in the rhizosphere soils of M+ plants (A, nodes = 296, edges =

853

4620, transitivity = 0.63) and M- plants (B, nodes = 309, edges = 3595, transitivity = 0.62),

854

respectively. A connection stands for a significant correlation (Spearman’s correlation test, PFDR <

855

0.05). Nodes belonging to Rhizobiales, Sphingomonadales, Actinomycetales, Xanthomonadales and

856

Burkholderiales members are colored. The size of each node is proportional to the degree (number of

857

connections) of each node.

858

(C) The distribution of degree centrality of bacterial orders across rhizosphere networks of M+ and

859

M- plants. Statistical significance was determined via multiple comparison with Kruskal-Wallis test

860

and with false discovery rate correction (α = 0.05).

861

(D and E) Co-abundance network of Spearman’s rho correlation between paired ASVs based on the

862

QMP data set of 16S rRNA genes in the roots of M+ plants (D, nodes = 264, edges = 3250,

863

transitivity = 0.46) and M- plants (E, nodes = 300, edges = 4705, transitivity = 0.59), respectively. A

864

connection stands for a significant correlation (Spearman’s correlation test, PFDR < 0.05). Nodes are

865

colored and sized as in A and B.

Jo

ur

na

lP

re

-p

ro

of

839

(F) The distribution of degree centrality of bacterial orders across root networks of M+ and M- plants.

867

Statistical significance was determined via multiple comparison with Kruskal-Wallis test and with

868

false discovery rate correction (α= 0.05).

869

Figure 4. AM Symbiosis Promotes Rhizobial Accumulation in the Rhizosphere of M. truncatula

870

and Efficient Rhizobium-Legume Nodulation in Native Soil.

871

(A) The load of rhizobia rpoB sequences across rhizosphere of plant genotypes. Data are mean ±

872

s.e.m., n = 5 biological replicates.

873

(B) Heatmap and phylogenetic tree of 14 rhizobia species retrieved by rpoB reads in rhizospheres of

874

different plant genotypes. The ML phylogenetic tree of rhizobia species was constructed using

875

representative rpoB gene sequences retrieved in this study (1000 bootstrap replications).

876

(C) Numbers of functional red nodules on roots of two-month-old plants grown in native soil. Data

877

are mean ± s.e.m., n = 12 biological replicates.

878

(D) Numbers of functional red nodules on roots of WT and pt4 plants grown in sterilized vermiculite

879

conditions (inoculated with rhizobia S. meliloti Rm2011 or S. medicae RS90, OD600 = 0.001). S.

880

medicae RS90, which share 100% rpoB gene sequence identify to the nodule isolates and pASV2

881

(Supplementary Figure 13), was isolated in rhizosphere soil of WT in this project. Data are mean ±

882

s.e.m., n = 12 biological replicates.

883

(E) Percentage of the nodulated P. vulgaris, P. sativum and G. soja plants inoculated with 1ml 10-5

884

dilutions of bulk soil, WT, lyk3, dmi3 and pt4 rhizosphere soil samples. n = 18 biological replicates.

885

For A-E, statistical significance was determined via multiple comparison with Kruskal-Wallis test

886

and with false discovery rate correction (α= 0.05). Post hoc test is indicated by letters at the top,

887

genotypes with the same letter are indistinguishable at 95% confidence.

888

Jo

ur

na

lP

re

-p

ro

of

866

*

Rhizosphere

M+R+ M+R-

0

Root

M-R- M-R+ M+R+ M+R-

M-R- M-R+

M+R+ M+R-

Rhizosphere

M-R- M-R+ M+R+ M+R-

* *

**

# # #

***

***

***

0.0e+00

*

*

* *

# 6.6e+8

2.2e+9

2.0e+9

1.1e+9

9.7e+8

9.2e+8

1.2e+9

** * 6.3e+9

4.4e+9

5.3e+9

8.0e+9

7.8e+9

#

2.5e+09 # #

5.0e+09

#

7.5e+09 8.9e+9

*

E

#

of

ro

-p

re

M-R+

# #

M-R+

lyk3_1 lyk3_3 lyk3_4 lyk3_2 lyk3_5 nfp_5 nfp_1 nfp_2 nfp_3 nfp_4 WT_4 WT_2 WT_5 WT_1 WT_3 ipd3_1 ipd3_2 ipd3_3 ipd3_4 ipd3_5 dmi2_5 dmi2_3 dmi2_1 pt4_3 pt4_2 pt4_4 pt4_1 pt4_5 dmi3_3 dmi3_5 dmi3_4 dmi3_1 dmi3_2 dmi2_4 dmi2_2

nfp ipd3 dmi2 dmi3 pt4

cluster4

0

#

M-Re

cluster5

cluster3

5.0e+5

6.4e+9

cluster6

M+R-

#

D

lP

cluster7

1.0e+6

#

ipd3 dmi2 dmi3 pt4

#

C Rhizobia only AMF only Adtively Non−additively

cluster2 cluster1

1.5e+6

QMP phyla abundance (g-1)

#

100

na

ur

0.05 nfp f

* * *

c

#

#

2.0e+6

#

#

0.10

WT lyk3 d

#

0.15 b

#

M-R-

# #

0.20 a a

c

#

* *

M+R-

c

#

#

50

Jo

ab

* * * * ***

M+R+ ab

# #

M+R+

#

75

* * *

a

# #

c

#

WT lyk3

#

#

Quantified AMF 18S abundance (g-1)

b

# #

#

#

Total phosphorus (g/100g)

B

# # #

*

25 #

* *

RMP 16S frequency (%)

A

scale log(FPKM） 2

1

0

−1

−2

Bacterial phyla

Alphaproteobacteria Actinobacteria Bacteroidetes Gammaproteobacteria Betaproteobacteria Acidobacteria Verrucomicrobia Chloroflexi Firmicutes Gemmatimonadetes Deltaproteobacteria Others

Root

M-R- M-R+

A

B Root microbiota ~ Genotype (33.2 % of variance; P = 1e-04) Constrained PCoA 2 (11.71%)

Constrained PCoA 2 (23.94%)

0.4

Rhizosphere microbiota ~ Genotype (43.7 % of variance; P = 1e-04)

0.0

-0.4

Genotype WT

0.4

lyk3 nfp ipd3

0.0

dmi2 dmi3 pt4

-0.4

-0.8 0.25

0.50

0.75

-0.2

Constrained PCoA 1 (69.45%)

*'+1&)$*'(%+$ 2+%/+* 3'(4

lP

65C86'%("'&/9*('$*($( 65CN8U( 65=7?8%%'%#$%*&&**34 )$*'(% 65 65 E=8%'#%&)$ 657B<88%'#%&)$*'(*'(% % 6 J!! 6 5?B8M +1 $*($(
na

;? 1 65 5==< E8++=P=3'34 '(% 6 5=< 8W"+' !+)$* 6 5B; 8A3'(& 6 5=N7 =N= %$*($( 6 5EE=8G &1+&)$*'( 6 >N8F#'# 65 E?B8J!!+1?=; 65 65778,%$-9%#+:&)+3/ 65;<8,%$65=>8,%$-99%%##++::&&)+3/ )+3/ 65? 65 ;8@#+:&)+$!( 65=?>=78A$)%94 4 6 8@# 6 5=?>8 &-&"!$ 65C5
'(% )$* % $3!& )$*'( E;8F $3!& 4+$ 65=5?78F?=8G(H&H&4+$ 6 5E C8G( 4+$ 6 = 65E5=E8G(H(&H&4+$ / 6 5E>8G %#+:&)+3 6 +1& / 65C==<868I(4&%#+:&)+3+3/ 65 (4&%#+:&) 65BE8II( +3/ 65?<8 4&%#+:&) 657<8J14+K(%

65 65 ?>8M 65> ;=8M$+4$+4'&) 65 C8M$ '&) $*'( 65==;N8M$+4'+&4'&)$*$'(*'(%% B8 ) 657;8M$+4'&)$*$*'(%% 65E8@@##+:&)+3/'(% 65?8 +:&)+3/ 65=?8@#+:@&#)+:+$&)+3/ 65EB8L%&)$*'(*%+(3$/( 65=78@#+:&)+$*($(

Jo

ur

65= 65 N<8I =N; ( 8G 65;'#9!&"#&S-&1 65 ?85$%+&+!$*($(!$ 657C=7<86FPQH&%$R( E P 65CC8F 65=7>8686FPQP>7>7 65=NE8F##%94(&)$*'(F%+P3QP>7 +'+1&"#$ / 65C?<8O+L$$4*'((!$( 65=E>8@#&-&4"+%+!!$*($!($ 65
!"#$"%&'(&)$*'(%+$

*+

Bacterial phyla Alphaproteobacteria Actinobacteria Bacteroidetes Gammaproteobacteria Betaproteobacteria Acidobacteria Verrucomicrobia Firmicutes Gemmatimonadetes

/+*%&)+$ %%3*& %+$5( $*'( -&)

re

-p

ro

+$ % ( ' * $ &) ( ' & % $" / .$/

65786'%("'&/9*(4 65>8@#&-&*7&8*F*=3=4= 65EN &1$4 %/&/ 65BC8UU##((%/&/&$1*$'4(% 65=;<857=8A94&&)/&1$4$4 6 #(%/ /&1 ( 8U & $ 65C7 8U#(%/1$-$*(*($( 4 65>E1'#&/& &1$-$&/&1$ $ / # N?8T 1'#& '%&" 65E =;N8T$ 86'(1& 65 65EE

,$*'(%&+,('$" %&'(& 0

C

0.2

of

Constrained PCoA 1 (50.98%)

0.0

.(/ /$'0

0.00

lyk3/WT nfp/WT ipd3/WT dmi2/WT dmi3/WT pt4/WT

-0.25

(-)log10FDR 2 1.3 0 -1.3 -2

A

D

M+ rhizosphere co−abundance network

M+ root co−abundance network

Bacterial orders Rhizobiales Sphingomonadales Actinomycetales Xanthomonadales Burkholderiales

B

of

Others

E

M- root co−abundance network

Jo

ur

na

lP

re

-p

ro

M- rhizosphere co−abundance network

C

1 00 <

P

= P

P

P

0.

70 0.

0. =

0. =

0. = P

4

7 11

6 00

0 92

4 0. = P

● ● ● ● ● ●

●

ia

le s

le s

●

● ● ● ● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ●

s

● ● ● ● ● ● ● ●

ad a

ad al

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

● ●

●

ho m Xa nt

on

●

●

● ● ●

● ● ● ● ● ●

● ● ● ●

th er

● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ●

O

●

on

s al e go m hi n

● ● ● ● ● ●

● ● ●

es

● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

● ● ●

● ● ● ● ● ● ● ●

er

● ● ● ● ● ● ● ● ● ●

●

● ●

ho ld

● ● ● ● ● ● ● ● ● ● ● ● ● ●

● ●

● ●

rk

● ● ●

●

Bu

● ●

● ● ● ● ● ●

bi Sp

th O

ia ld er

rk ho

0

er s

le s

●

● ● ● ● ● ●

● ● ● ● ● ● ● ●

● ● ● ●

hi zo

● ● ● ● ●

R

● ● ● ● ●

● ● ●

● ● ● ● ● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ●

● ●

●

● ●

● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

20

●

es

● ● ●

● ● ● ● ● ● ● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

et al

● ● ● ● ● ● ● ●

40

●

● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

Symbiosis

●

● ●

om yc

●

● ●

● ● ● ●

60

tin

●

10

4 =

P ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

Degree centrality of root networks

●

● ● ●

● ●

● ● ● ● ● ●

● ●

Ac

● ● ● ● ● ● ● ● ● ● ● ● ● ●

● ●

● ● ●

●

●

● ● ● ● ● ● ● ●

● ●

0.

42 0. = P

● ● ● ● ● ●

● ● ●

● ● ● ● ●

01

8

1 0. P

●

● ●

●

ho m nt Xa

hi Sp

● ●

● ●

al e

m

R

tin

02

9 ● ● ● ● ●

on ad

zo bi hi

●

● ● ●

le s

● ● ● ● ● ●

al

al es

● ● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

es

● ●

● ●

● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ●

● ● ●

● ● ● ● ● ●

● ●

● ● ●

s

●

● ● ● ● ● ● ● ● ● ●

● ● ●

●

on a

● ●

●

● ● ●

Bu

● ● ● ● ● ● ● ● ●

● ●

● ● ● ● ● ● ● ● ● ● ●

●

da

● ● ●

● ● ● ● ● ●

=

0. ● ● ● ● ●

● ●

0

Ac

P

●

●

● ● ● ● ●

● ● ● ● ● ●

ng o

20

●

● ● ●

● ● ● ●

●

40

=

0. = P ● ● ● ● ● ● ●

● ●

● ●

60

06

9 00

1 48 0. = P 80

om yc et

Degree centrality of rhizosphere networks

F

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

●

M+ plants

●

M− plants

Red nodule numbers per plant in sterile sand 5

4

3

2

-p

a

a b c

c

a

3e+08

c

c b

c

0

Phaseolus vulgaris

b

c

b

b

5e+08

M+R+ M+RM-R-

c b

a

D

6

4

a

c

10.0

bc

10

0

Pisum sativum

of

a R. anhuiense R. mesosinicum P. giardinii S. medicae M. plurifarium M. huakuii M. silamurunense B. elkanii B. jicamae B. lablabi B. yuanmingense B. huanghuaihaiense B. diazoefficiens B. japonicum

B

ro

a

Red nodule numbers per plant in sterile vermiculite

a

na

2

re

0

lP

E

ur

Quantified rhizobia rpoB abundance (g-1)

C

Jo

Red nodule numbers per plant in native soil

A

C

1e+08

a

c

a

a

a

7.5

5.0

2.5

0.0

b c bc

6

20

4

1 2

0

Glycine soja c

(-)log10FDR 2 1 0 -1 -2

lyk3/WT nfp/WT ipd3/WT dmi2/WT dmi3/WT pt4/WT

5e+07

a Strain Sm1021 RS90