The rhizospheric bacterial diversity of 
Chuan BeiMu is a kind of the precious traditional Chinese herbal medicine, its source materials are the bulbs of 5 species in genus Fritillaria ( Fritillaria cirrhosa D. Don, Fritillaria przewalski i Maxim., Fritillaria unibracteata Hsiao et K. C. Hsia, Fritillaria delavayi Franch. and Fritillaria taipaiensis P. Y. Li) . The Chuan Beimu medicine contain multiple compounds such as alkaloids, steroids, terpenoids, purine and fatty acid, which can cool down body temperature, relieve coughing and asthma, and resolve phlegm – . Fritillaria taipaiensis P. Y. Li, which has been collected in Chinese Pharmacopoeia as one of the medical material provenances of Chuan BeiMu in 2015 , has important economic value. In recent years, the production of F. taipaiensis has been popularized in suitable areas – . However, in the actual production, the seedlings of F. taipaiensis may appear the problems such as growth potential reducing , disease accumulation , and quality decline after 3–4 years’ continuous planting, which badly affected its production and quality. In response to this situation, current researches mainly focus on the fertilization control and the planting technology improvement , . Soil contains different and wide variety of microbial communities, and the interaction between soil microorganism and soil factors plays an important role in the process of soil organic matter degradation and soil nutrient cycling , . Meanwhile, the quality of soil can also be reflected by the soil enzyme activities and the microbial community composition , . Researches have demonstrated continuous cropping and long-term single plating pattern could affect the diversity and structure of the microorganism population in the planting soil – , these microorganism and then could directly or indirectly affect the growth of plants by changing the chemical properties of soil . In addition, the pathogenic microorganisms accumulated during the long-term single planting can also hinder the growth of plants – . The studies on the changes of rhizospheric microbial community structure during the F. taipaiensis cultivation have not been deeply progressed at present. In this study, the population structure compositions, richness and evenness of the bacteria in the rhizospheric soils, in which F. taipaiensis was planted for 1–5 years, were compared and analyzed to reveal the change on the bacterial diversity in the rhizosphere of Fritillaria taipaiensis under long-term single planting. Moreover, we analyzed the relationship between the soil bio-chemical factors and the rhizospheric bacterial community to explore the interaction between soil environment and soil microorganisms. The result of this study could provide a theoretical basis for the soil management in the sustainable cultivation of F. taipaiensis . 

Sampling site location The sampling site is a private planting base located in Lanying Village, Wuxi County, Chongqing, China (31°23′56.11″ N, 109°50′29.93″ E-31°35′26.57″ N, 109°00′11.96″ E; 2 274–2 290 m above sea level), with the annual average temperature of 7.2 °C, and the annual average precipitation of 1 100–1 300 mm, the soil in the study site was classified as mountain yellow–brown soil (fine-loamy, mixed, mesic Aridic Haplustalf) . The field was 20 × 20 m for the F. taipaiensis with each cultivation year, The plant and row spacing of F. taipaiensis was 15 × 15 cm, and the fields were uncultivated wastelands with the same soil basement condition before F. taipaiensis planting. The fertilization of the field was N-P-K 15:15:15, 225 kg hm −2 for each year. Sample collection In May, 2017, with the permission of the planting base owner, F. taipaiensis rhizospheric soils with different cultivation years from 1 to 5 (Y1: rhizospheric soil with 1 year of cultivation, Y2: rhizospheric soil with 2 years of cultivation, Y3: rhizospheric soil with 3 years of cultivation, Y4: rhizospheric soil with 4 years of cultivation, Y5: rhizospheric soil with 5 years of cultivation) were sampled. 3 sampling plots (3 × 3 m) were randomly selected in the field with each cultivation year. The regolith, impurities and gravels on the surface of fields were removed to expose the roots, and the soils attached on the roots 6 ~ 10 cm below the ground were collected by shaking root method . Soils from the rhizosphere of 5 F. taipaiensis plants randomly selected in each sampling plot were collected and then were mixed. Then the soil samples were placed into the aseptic bags in ice box and then were transferred into laboratory after removed the litters and residual cover. The experimental samples were sifted by 2 mm sieve and then were stored in − 80 °C ultra-low temperature freezer after freeze drying to be used in soil total DNA extraction. The bio-chemical properties of the F. taipaiensis rhizospheric soils with different cultivation years were shown in Table . Soil DNA extraction and high throughput sequencing 0.25 g F. taipaiensis rhizospheric soil sample was used for DNA extraction by CTAB method . The genomic DNA after diluted to 1 ng μl -1 was amplified by PCR via high-fidelity DNA polymerase and Phusion® High-Fidelity PCR Master Mix with GC Buffer (New England Biolabs) by using specific primer 515F (5′-GTGCCACCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) for the bacterial 16sRNA gene V4 sequence, the amplification condition was pre-denaturation at 98 °C for 2 min, denaturation at 98 °C for 30 s, annealing at 50 °C for 30 s, extending at 72 °C for 1 min, after 30 cycles, and then extending 72 °C again for 5 min to end the PCR reaction. The PCR products were tested by 2% agarose gel electrophoresis and then were extracted by gel extraction kit (Qiagen). The DNA library was constructed by using TruSeq® DNA PCR-Free Sample Preparation Kit (Illumina) and then was sequenced by HiSeq2500 PE250. Data analysis The data of each sample was separated from the off machine data according to the Barcode sequence and PCR amplification primer sequence, and the reads of each sample were assembled by FLASH (V1.2.7, http://ccb.jhu.edu/software/FLASH/ ) after removed the Barcode sequence and the PCR amplification primer sequence to obtain the Raw date. The Raw date were filtered and spliced to get Clean date refer to the Tags quality control process in QIIME (V1.9.1, http://qiime.org/scripts/split_libraries_fastq.html ). And then, the chimeras in the Clean date were eliminated by comparing with UCHIME (Algorithm, http://www.drive5.com/usearch/manual/uchime_algo.html ) and database (Unite database, https://unite.ut.ee/ ) to obtain the Effective tags. The depth of sequencing was 40 000 reads at least in each original library, the operational taxonomic unit (OTU) was analyzed and clustered by using the UPARSE software (version 7.0.1090) based on 97% similarity. According to the result of OTU clustering, the representative sequence of each OTU was annotated by a species by using QIIME (version 1.9.1). The relative abundance was analyzed via Mothur (version.1.44.0) in the SS rRNA database of SILVA ( http://arb-silva.de/ ). The Mothur software (version.1.44.0) was used to calculate the Goods-coverage, Chao1, ACE, Shannon and Simpson indices. The QIIME software (version 1.9.1) was used to calculate the Unifrac distance and plot the unweighted pair-group method using arithmetic average (UPGMA) cluster tree. The R-software (version 2.15.3) was used to plot the Rank abundance curve. Meanwhile, the relative abundance, the Alpha diversity of OTUs was analyzed to get the information such as the unique species abundance in one sample, the common and the unique OTUs among the different samples and groups by Tukey-test. Non-Metric Multi-Dimensional Scaling (NMDS) analysis, distance-based redundancy analysis (db-RDA) and permutational MANOVA analysis (Adonis) were performed by vegan bag on R-software. The data were analyzed by single factor (one-way ANOVA) and Duncan for variance analysis and multiple comparison via SPSS 22.0 and Excel 2003 (α = 0.05). Correlation analysis was conducted on SPSS 22.0 using the Pearson method (α = 0.05). 

The sampling site is a private planting base located in Lanying Village, Wuxi County, Chongqing, China (31°23′56.11″ N, 109°50′29.93″ E-31°35′26.57″ N, 109°00′11.96″ E; 2 274–2 290 m above sea level), with the annual average temperature of 7.2 °C, and the annual average precipitation of 1 100–1 300 mm, the soil in the study site was classified as mountain yellow–brown soil (fine-loamy, mixed, mesic Aridic Haplustalf) . The field was 20 × 20 m for the F. taipaiensis with each cultivation year, The plant and row spacing of F. taipaiensis was 15 × 15 cm, and the fields were uncultivated wastelands with the same soil basement condition before F. taipaiensis planting. The fertilization of the field was N-P-K 15:15:15, 225 kg hm −2 for each year. 

In May, 2017, with the permission of the planting base owner, F. taipaiensis rhizospheric soils with different cultivation years from 1 to 5 (Y1: rhizospheric soil with 1 year of cultivation, Y2: rhizospheric soil with 2 years of cultivation, Y3: rhizospheric soil with 3 years of cultivation, Y4: rhizospheric soil with 4 years of cultivation, Y5: rhizospheric soil with 5 years of cultivation) were sampled. 3 sampling plots (3 × 3 m) were randomly selected in the field with each cultivation year. The regolith, impurities and gravels on the surface of fields were removed to expose the roots, and the soils attached on the roots 6 ~ 10 cm below the ground were collected by shaking root method . Soils from the rhizosphere of 5 F. taipaiensis plants randomly selected in each sampling plot were collected and then were mixed. Then the soil samples were placed into the aseptic bags in ice box and then were transferred into laboratory after removed the litters and residual cover. The experimental samples were sifted by 2 mm sieve and then were stored in − 80 °C ultra-low temperature freezer after freeze drying to be used in soil total DNA extraction. The bio-chemical properties of the F. taipaiensis rhizospheric soils with different cultivation years were shown in Table . 

0.25 g F. taipaiensis rhizospheric soil sample was used for DNA extraction by CTAB method . The genomic DNA after diluted to 1 ng μl -1 was amplified by PCR via high-fidelity DNA polymerase and Phusion® High-Fidelity PCR Master Mix with GC Buffer (New England Biolabs) by using specific primer 515F (5′-GTGCCACCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) for the bacterial 16sRNA gene V4 sequence, the amplification condition was pre-denaturation at 98 °C for 2 min, denaturation at 98 °C for 30 s, annealing at 50 °C for 30 s, extending at 72 °C for 1 min, after 30 cycles, and then extending 72 °C again for 5 min to end the PCR reaction. The PCR products were tested by 2% agarose gel electrophoresis and then were extracted by gel extraction kit (Qiagen). The DNA library was constructed by using TruSeq® DNA PCR-Free Sample Preparation Kit (Illumina) and then was sequenced by HiSeq2500 PE250. 

The data of each sample was separated from the off machine data according to the Barcode sequence and PCR amplification primer sequence, and the reads of each sample were assembled by FLASH (V1.2.7, http://ccb.jhu.edu/software/FLASH/ ) after removed the Barcode sequence and the PCR amplification primer sequence to obtain the Raw date. The Raw date were filtered and spliced to get Clean date refer to the Tags quality control process in QIIME (V1.9.1, http://qiime.org/scripts/split_libraries_fastq.html ). And then, the chimeras in the Clean date were eliminated by comparing with UCHIME (Algorithm, http://www.drive5.com/usearch/manual/uchime_algo.html ) and database (Unite database, https://unite.ut.ee/ ) to obtain the Effective tags. The depth of sequencing was 40 000 reads at least in each original library, the operational taxonomic unit (OTU) was analyzed and clustered by using the UPARSE software (version 7.0.1090) based on 97% similarity. According to the result of OTU clustering, the representative sequence of each OTU was annotated by a species by using QIIME (version 1.9.1). The relative abundance was analyzed via Mothur (version.1.44.0) in the SS rRNA database of SILVA ( http://arb-silva.de/ ). The Mothur software (version.1.44.0) was used to calculate the Goods-coverage, Chao1, ACE, Shannon and Simpson indices. The QIIME software (version 1.9.1) was used to calculate the Unifrac distance and plot the unweighted pair-group method using arithmetic average (UPGMA) cluster tree. The R-software (version 2.15.3) was used to plot the Rank abundance curve. Meanwhile, the relative abundance, the Alpha diversity of OTUs was analyzed to get the information such as the unique species abundance in one sample, the common and the unique OTUs among the different samples and groups by Tukey-test. Non-Metric Multi-Dimensional Scaling (NMDS) analysis, distance-based redundancy analysis (db-RDA) and permutational MANOVA analysis (Adonis) were performed by vegan bag on R-software. The data were analyzed by single factor (one-way ANOVA) and Duncan for variance analysis and multiple comparison via SPSS 22.0 and Excel 2003 (α = 0.05). Correlation analysis was conducted on SPSS 22.0 using the Pearson method (α = 0.05). 

Alpha diversity analysis The Rank abundance curves (Fig. ) which were drawn according to the relative abundance of the OTUs showed: the curve of sample Y2 had the maximum abscissa span, which indicated that the F. taipaiensis rhizospheric soil with 2 years of cultivation had the maximum bacterial species richness; and the curve of sample Y1 had the softest changing trend in ordinate which indicated the F. taipaiensis rhizospheric soil with 1 year of cultivation had the maximum bacterial species evenness; the curve of sample Y4 had the minimum abscissa span and the maxim ordinate changing trend which indicated the F. taipaiensis rhizospheric soil with 4 years of cultivation had the minimum bacterial species richness and evenness. As shown in Table : Goods_coverage showed the species coverage detected in each sample was above 97%; there was no significant difference in the Chao1 and ACE indices between the samples, which reflected the species richness of the samples with different cultivation years had no significant change; meanwhile, the Shannon and Simpson indices of the samples showed a trend of declining with the cultivation year, which indicated the reduction of bacterial community diversity. Beta diversity analysis of the bacterial communities in the rhizospheric soils As shown in Fig. , the bacterial phyla with the top 10 relative abundances in the F. taipaiensis rhizospheric soils with different cultivation years were Proteobacteria (34.16–54.82%), Acidobacteria (10.19–18.16%), Bacteroidetes (9.49–13.75%), Verrucomicrobia (4.91%-7.21%), Firmicutes (0.18–6.05%), Gemmatimonadetes (2.87–5.93%), Actinobacteria (1.51–4.55%), Planctomycetes (2.05–3.06%), Nitrospirae (1.01–2.48%) and Chloroflexi (0.96–1.83%). The dominant phylum was Proteobacteria which presented a rising relative abundance from Y1 to Y4 with the cultivation year. And the relative abundance of Firmicutes, Actinobacteria and Chloroflexi showed a downward changing trend from Y1 to Y4, in which the relative abundance of Firmicutes had great variation among samples, declined from 6.05% in Y1 to 0.65% in Y3. Meanwhile, the samples were clustered by unweighted pair group method with arithmetic means (UPGMA). Samples Y4 and Y5 could be grouped together but be separated from the cluster of samples Y1, Y2 and Y3, which revealed that the samples with 4 and 5 years of cultivation had the similar bacterial community structures but were different from the samples with 1, 2, and 3 years of cultivation. The samples Y1, Y2, and Y3 could be grouped in a large cluster, in which the samples with 1 and 3 years of cultivation could be further grouped in a smaller cluster. The bacterial genera with the top 15 relative abundances in the F. taipaiensis rhizospheric soil were Methylotenera , Pseudomonas , Lactobacillus , Gemmatimonas , Ralstonia , Bryobacter , Sphingobium , Flavobacterium , Collimonas , Sphingomonas , Candidatus , Rhizomicrobium , Pedobacter , Bacteroides and unculturable bacterium RB41 (Table ). The relative abundances of Methylotenera , Pseudomonas and Sphingobium showed a tendency of ascending with the cultivation year, the relative abundance of Methylotenera in Y4 was 16.07 times of that in Y1, and the relative abundance of Pseudomonas and Sphingobium in Y5 was 5.21 and 4.15 times of that in Y1, respectively. The relative abundances of Lactobacillus , Gemmatimonas , Bryobacter and Bacteroides declined with the cultivation year, the relative abundance of Gemmatimonas , Bryobacter and Bacteroides in Y5 declined by 4.67, 2.29, and 99 times compared with that in Y1, respectively, and the relative abundances of Lactobacillus in Y4 and Y5 were absent. The relative abundances of Flavobacterium , Collimonas and Pedobacter showed a downward trend after a initial rising, the relative abundance of Flavobacterium in Y4 and the relative abundance of Collimonas in Y3 was 5.89 and 12 times of that in Y1, respectively. While in Y5, the relative abundances of these 2 genera declined significantly when compared with those in Y3 and Y4, respectively. The relative abundances of Pedobacter had no obvious change in Y1, Y2 and Y3 but rose suddenly to a high level in Y4 and then declined obviously in Y5. Meanwhile, the relative abundances of Candidatus also had no significant change in Y1, Y2 and Y3 but declined obviously in Y4 and Y5. All the results above indicated that F. taipaiensis rhizospheric soils with different cultivation years had significant difference in the dominate bacteria community on the genus level. The heat map of the bacterial genera with the top 35 relative abundances in each sample showed a significant difference between the dominate bacterial communities in the F. taipaiensis rhizospheric soils with different cuntivation years (Fig. A). The dominant bacterial genera in F. taipaiensis rhizospheric soil with different cultivation years were Y1: Haliangium , Ferruginibacter , Chthoniobacter , Germmatimonas , Bacteroides , Lactobacillus , Opitutus , Arenimonas and Piscinibacter ; Y2: Pseudarthrobacter , Lachnospiraceae_UCG_001 and Sphingomonas ; Y3: Steroidobacter , Duganella , Collimonas , Aeromonas , unidentified Acidobacteriaceae_(Subgroup_1) and Mucilaginibacter ; Y4: H16, Methylotenera , Pedobacter , Lysobacter , Polycyclovorans , Massilia and Flavobacterium ; Y5: Duganella , Pseudomonas , Pseudarthrobacter , Haliangium , RB41 and Acidibacter . As shown in the Venn diagram, the samples with different cultivation years shared 2 928 common OTUs, the number of the unique OTUs in each sample from Y1 to Y5 was 339, 490, 213, 281 and 332, respectively (Fig. B). It indicated that the sample Y2 was greatly different with the other samples in bacterial community structures. Heat map was created via QIIME (V1.9.1, http://qiime.org/scripts/split_libraries_fastq.html ). The values corresponding to the heat map are Z-values obtained by standardizing the relative abundances of the species in each row. The Z-value of a sample in a certain classification means the difference between the relative abundance of the sample and the average relative abundance of all samples divided by the standard deviation of all samples in the classification. NMDS and adonis analysis As shown in Fig. , the points representing samples Y1 and Y3 overlapped in their distribution ranges to a certain degree. However, the points representing these two samples were far from the other samples, indicating the rhizospheric bacterial communities in Y1 and Y3 were similar but were significantly different from the others. The above result was similar to that of UPGMA clustering. The further analysis by Adonis showed there were significant differences in rhizospheric bacterial communities between Y1 and Y2, Y2 and Y3, Y1 and Y5, Y3 and Y5 ( p = 0.001) (Table ). Relationship between the rhizospheric bacterial community and the soil bio-chemical factors The correlation analysis of the soil bio-chemical factors showed: the pH was significantly negatively correlated with the contents of available phosphorus, available nitrogen, organic matter, and the activity of urease ( p < 0.01), was significantly positively correlated with the activity of invertase ( p < 0.01), and was negatively correlated with the activity of alkaline phosphatase ( p < 0.05); the activity of acid phosphatase was significantly positively correlated with the content of available potassium ( p < 0.01); the activity of alkaline phosphatase was significantly positively correlated with the contents of available phosphorus, available nitrogen, organic matter and the activity of urease ( p < 0.01), and was significantly negatively correlated with the activities of protease and invertase ( p < 0.01); the urease activity was significantly positively correlated with the contents of available phosphorus, available nitrogen and organic matter ( p < 0.01), was significantly negatively correlated with the activity of invertase ( p < 0.01), and was negatively correlated with the protease activity ( p < 0.05); the activity of invertase was significantly negatively correlated with the contents of available phosphorus, available nitrogen, organic matter, and was significantly positively correlated with the activity of protease ( p < 0.01); the activity of protease was significantly negatively correlated with the contents of available phosphorus and available nitrogen ( p < 0.01), and was negatively correlated with the content of organic matter ( p < 0.05); the activity of catalase was positively correlated with the content of organic matter ( p < 0.05); the content of organic matter was significantly positively correlated with the content of available phosphorus ( p < 0.01), and was positively correlated with the content of available nitrogen ( p < 0.05); the content of available nitrogen was significantly positively correlated with the content of available phosphorus ( p < 0.01) (Table ). Another correlation was established between the soil bio-chemical factors and the OTU numbers of the bacterial genera with the top 15 relative abundances: the pH was significantly positively correlated with the abundances of Methylotenera , Ralstonia and Flavobacterium ( p < 0.01), was positively correlated with the abundance of Pedobacter ( p < 0.05), was significantly negatively correlated with the abundances of Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), and was negatively correlated with the abundances of Lactobacillus and Gemmatimonas ( p < 0.05); the acid phosphatase activity was significantly negatively correlated with the abundance of Flavobacterium ( p < 0.01), negatively correlated with the abundances of Methylotenera and Collimonas ( p < 0.05), and was positively correlated with the abundances of Gemmatimonas and RB41 ( p < 0.05); the alkaline phosphatase activity was significantly positively correlated with the abundances of Lactobacillus , Gemmatimonas and Bryobacter ( p < 0.01), was significantly negatively correlated with the abundances of Pseudomonas and Sphingobium ( p < 0.01), was positively correlated with the abundances of Candidatus , Rhizomicrobium and Bacteroides ( p < 0.05), and was negatively correlated with the abundances of Methylotenera and Ralstonia ( p < 0.05); the urease activity was significantly positively correlated with the abundances of Lactobacillus , Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), was positively correlated with the abundance of Gemmatimonas and Bacteroides ( p < 0.05), was significantly negatively with the abundances of Methylotenera and Flavobacterium ( p < 0.01), and was negatively correlated with the abundance of Ralstonia ( p < 0.05); the invertase was significantly positively correlated with the abundances of Ralstonia and RB41 ( p < 0.01), was positively correlated with the abundances of Methylotenera and Pseudomonas ( p < 0.05), was significantly negatively correlated with the abundances of Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), and was negatively correlated with the abundance of Lactobacillus ( p < 0.05); the protease activity was significantly positively correlated with the abundances of Pseudomonas and RB41 ( p < 0.01), was positively correlated with the abundance of Ralstonia and Collimonas ( p < 0.05), was significantly negatively correlated with the abundance of Bryobacter ( p < 0.01), and was negatively correlated with the abundances of Lactobacillus , Candidatus and Rhizomicrobium ( p < 0.05); the catalase activity was positively correlated with the abundance of Gemmatimonas and was negatively correlated with the abundance of Sphingomonas ( p < 0.05); the available nitrogen content was significantly positively correlated with the abundances of Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), was positively correlated with the abundance of Lactobacillus and Sphingomonas ( p < 0.05), and was negatively correlated with the abundances of Methylotenera , Ralstonia and Flavobacterium ( p < 0.05); the available phosphorus content was significantly positively correlated with the abundances of Lactobacillus , Gemmatimonas , Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), was positively correlated with the abundance of Bacteroides ( p < 0.05), was significantly negatively correlated with the abundances of Methylotenera , Ralstonia and Sphingobium ( p < 0.01), and was negatively correlated with the abundances of Pseudomonas and Flavobacterium ( p < 0.05); the content of organic matter was significantly positively correlated with the abundances of Lactobacillus and Gemmatimonas ( p < 0.01), was positively correlated with the abundances of Bryobacter , Candidatus and Rhizomicrobium ( p < 0.05), was significantly negatively correlated with the abundances of Methylotenera , Pseudomonas , Ralstonia and Sphingobium ( p < 0.01), and was negatively correlated with the abundance of Flavobacterium ( p < 0.05) (Table ). The results of the db-RDA showed the effect of the bacterial community on the soil bio-chemical factors in the F. taipaiensis rhizospheric soils with different cultivation years explained 68.36% of the total characteristic value (Fig. A). Further analysis showed that the pH (R 2 = 0.9005, p = 0.001), the urease activity (R 2 = 0.7792, p = 0.001), the available phosphorus content (R 2 = 0.7890, p = 0.001) and the organic matter content (R 2 = 0.8564, p = 0.001) were the dominate factors that affected by the change of the rhizospheric soil bacterial community structure, among which the pH was the most significant factor (Supplementary Table ). The average distance from the point to the center in each group of the PCoA ranking result based on the Bray–curtis distance matrix could represent the beta diversity within group. And the beta diversity of the bacterial community within group showed an upward trend after an initial decrease with the cultivation year (Fig. B). 

The Rank abundance curves (Fig. ) which were drawn according to the relative abundance of the OTUs showed: the curve of sample Y2 had the maximum abscissa span, which indicated that the F. taipaiensis rhizospheric soil with 2 years of cultivation had the maximum bacterial species richness; and the curve of sample Y1 had the softest changing trend in ordinate which indicated the F. taipaiensis rhizospheric soil with 1 year of cultivation had the maximum bacterial species evenness; the curve of sample Y4 had the minimum abscissa span and the maxim ordinate changing trend which indicated the F. taipaiensis rhizospheric soil with 4 years of cultivation had the minimum bacterial species richness and evenness. As shown in Table : Goods_coverage showed the species coverage detected in each sample was above 97%; there was no significant difference in the Chao1 and ACE indices between the samples, which reflected the species richness of the samples with different cultivation years had no significant change; meanwhile, the Shannon and Simpson indices of the samples showed a trend of declining with the cultivation year, which indicated the reduction of bacterial community diversity. 

As shown in Fig. , the bacterial phyla with the top 10 relative abundances in the F. taipaiensis rhizospheric soils with different cultivation years were Proteobacteria (34.16–54.82%), Acidobacteria (10.19–18.16%), Bacteroidetes (9.49–13.75%), Verrucomicrobia (4.91%-7.21%), Firmicutes (0.18–6.05%), Gemmatimonadetes (2.87–5.93%), Actinobacteria (1.51–4.55%), Planctomycetes (2.05–3.06%), Nitrospirae (1.01–2.48%) and Chloroflexi (0.96–1.83%). The dominant phylum was Proteobacteria which presented a rising relative abundance from Y1 to Y4 with the cultivation year. And the relative abundance of Firmicutes, Actinobacteria and Chloroflexi showed a downward changing trend from Y1 to Y4, in which the relative abundance of Firmicutes had great variation among samples, declined from 6.05% in Y1 to 0.65% in Y3. Meanwhile, the samples were clustered by unweighted pair group method with arithmetic means (UPGMA). Samples Y4 and Y5 could be grouped together but be separated from the cluster of samples Y1, Y2 and Y3, which revealed that the samples with 4 and 5 years of cultivation had the similar bacterial community structures but were different from the samples with 1, 2, and 3 years of cultivation. The samples Y1, Y2, and Y3 could be grouped in a large cluster, in which the samples with 1 and 3 years of cultivation could be further grouped in a smaller cluster. The bacterial genera with the top 15 relative abundances in the F. taipaiensis rhizospheric soil were Methylotenera , Pseudomonas , Lactobacillus , Gemmatimonas , Ralstonia , Bryobacter , Sphingobium , Flavobacterium , Collimonas , Sphingomonas , Candidatus , Rhizomicrobium , Pedobacter , Bacteroides and unculturable bacterium RB41 (Table ). The relative abundances of Methylotenera , Pseudomonas and Sphingobium showed a tendency of ascending with the cultivation year, the relative abundance of Methylotenera in Y4 was 16.07 times of that in Y1, and the relative abundance of Pseudomonas and Sphingobium in Y5 was 5.21 and 4.15 times of that in Y1, respectively. The relative abundances of Lactobacillus , Gemmatimonas , Bryobacter and Bacteroides declined with the cultivation year, the relative abundance of Gemmatimonas , Bryobacter and Bacteroides in Y5 declined by 4.67, 2.29, and 99 times compared with that in Y1, respectively, and the relative abundances of Lactobacillus in Y4 and Y5 were absent. The relative abundances of Flavobacterium , Collimonas and Pedobacter showed a downward trend after a initial rising, the relative abundance of Flavobacterium in Y4 and the relative abundance of Collimonas in Y3 was 5.89 and 12 times of that in Y1, respectively. While in Y5, the relative abundances of these 2 genera declined significantly when compared with those in Y3 and Y4, respectively. The relative abundances of Pedobacter had no obvious change in Y1, Y2 and Y3 but rose suddenly to a high level in Y4 and then declined obviously in Y5. Meanwhile, the relative abundances of Candidatus also had no significant change in Y1, Y2 and Y3 but declined obviously in Y4 and Y5. All the results above indicated that F. taipaiensis rhizospheric soils with different cultivation years had significant difference in the dominate bacteria community on the genus level. The heat map of the bacterial genera with the top 35 relative abundances in each sample showed a significant difference between the dominate bacterial communities in the F. taipaiensis rhizospheric soils with different cuntivation years (Fig. A). The dominant bacterial genera in F. taipaiensis rhizospheric soil with different cultivation years were Y1: Haliangium , Ferruginibacter , Chthoniobacter , Germmatimonas , Bacteroides , Lactobacillus , Opitutus , Arenimonas and Piscinibacter ; Y2: Pseudarthrobacter , Lachnospiraceae_UCG_001 and Sphingomonas ; Y3: Steroidobacter , Duganella , Collimonas , Aeromonas , unidentified Acidobacteriaceae_(Subgroup_1) and Mucilaginibacter ; Y4: H16, Methylotenera , Pedobacter , Lysobacter , Polycyclovorans , Massilia and Flavobacterium ; Y5: Duganella , Pseudomonas , Pseudarthrobacter , Haliangium , RB41 and Acidibacter . As shown in the Venn diagram, the samples with different cultivation years shared 2 928 common OTUs, the number of the unique OTUs in each sample from Y1 to Y5 was 339, 490, 213, 281 and 332, respectively (Fig. B). It indicated that the sample Y2 was greatly different with the other samples in bacterial community structures. Heat map was created via QIIME (V1.9.1, http://qiime.org/scripts/split_libraries_fastq.html ). The values corresponding to the heat map are Z-values obtained by standardizing the relative abundances of the species in each row. The Z-value of a sample in a certain classification means the difference between the relative abundance of the sample and the average relative abundance of all samples divided by the standard deviation of all samples in the classification. 

As shown in Fig. , the points representing samples Y1 and Y3 overlapped in their distribution ranges to a certain degree. However, the points representing these two samples were far from the other samples, indicating the rhizospheric bacterial communities in Y1 and Y3 were similar but were significantly different from the others. The above result was similar to that of UPGMA clustering. The further analysis by Adonis showed there were significant differences in rhizospheric bacterial communities between Y1 and Y2, Y2 and Y3, Y1 and Y5, Y3 and Y5 ( p = 0.001) (Table ). 

The correlation analysis of the soil bio-chemical factors showed: the pH was significantly negatively correlated with the contents of available phosphorus, available nitrogen, organic matter, and the activity of urease ( p < 0.01), was significantly positively correlated with the activity of invertase ( p < 0.01), and was negatively correlated with the activity of alkaline phosphatase ( p < 0.05); the activity of acid phosphatase was significantly positively correlated with the content of available potassium ( p < 0.01); the activity of alkaline phosphatase was significantly positively correlated with the contents of available phosphorus, available nitrogen, organic matter and the activity of urease ( p < 0.01), and was significantly negatively correlated with the activities of protease and invertase ( p < 0.01); the urease activity was significantly positively correlated with the contents of available phosphorus, available nitrogen and organic matter ( p < 0.01), was significantly negatively correlated with the activity of invertase ( p < 0.01), and was negatively correlated with the protease activity ( p < 0.05); the activity of invertase was significantly negatively correlated with the contents of available phosphorus, available nitrogen, organic matter, and was significantly positively correlated with the activity of protease ( p < 0.01); the activity of protease was significantly negatively correlated with the contents of available phosphorus and available nitrogen ( p < 0.01), and was negatively correlated with the content of organic matter ( p < 0.05); the activity of catalase was positively correlated with the content of organic matter ( p < 0.05); the content of organic matter was significantly positively correlated with the content of available phosphorus ( p < 0.01), and was positively correlated with the content of available nitrogen ( p < 0.05); the content of available nitrogen was significantly positively correlated with the content of available phosphorus ( p < 0.01) (Table ). Another correlation was established between the soil bio-chemical factors and the OTU numbers of the bacterial genera with the top 15 relative abundances: the pH was significantly positively correlated with the abundances of Methylotenera , Ralstonia and Flavobacterium ( p < 0.01), was positively correlated with the abundance of Pedobacter ( p < 0.05), was significantly negatively correlated with the abundances of Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), and was negatively correlated with the abundances of Lactobacillus and Gemmatimonas ( p < 0.05); the acid phosphatase activity was significantly negatively correlated with the abundance of Flavobacterium ( p < 0.01), negatively correlated with the abundances of Methylotenera and Collimonas ( p < 0.05), and was positively correlated with the abundances of Gemmatimonas and RB41 ( p < 0.05); the alkaline phosphatase activity was significantly positively correlated with the abundances of Lactobacillus , Gemmatimonas and Bryobacter ( p < 0.01), was significantly negatively correlated with the abundances of Pseudomonas and Sphingobium ( p < 0.01), was positively correlated with the abundances of Candidatus , Rhizomicrobium and Bacteroides ( p < 0.05), and was negatively correlated with the abundances of Methylotenera and Ralstonia ( p < 0.05); the urease activity was significantly positively correlated with the abundances of Lactobacillus , Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), was positively correlated with the abundance of Gemmatimonas and Bacteroides ( p < 0.05), was significantly negatively with the abundances of Methylotenera and Flavobacterium ( p < 0.01), and was negatively correlated with the abundance of Ralstonia ( p < 0.05); the invertase was significantly positively correlated with the abundances of Ralstonia and RB41 ( p < 0.01), was positively correlated with the abundances of Methylotenera and Pseudomonas ( p < 0.05), was significantly negatively correlated with the abundances of Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), and was negatively correlated with the abundance of Lactobacillus ( p < 0.05); the protease activity was significantly positively correlated with the abundances of Pseudomonas and RB41 ( p < 0.01), was positively correlated with the abundance of Ralstonia and Collimonas ( p < 0.05), was significantly negatively correlated with the abundance of Bryobacter ( p < 0.01), and was negatively correlated with the abundances of Lactobacillus , Candidatus and Rhizomicrobium ( p < 0.05); the catalase activity was positively correlated with the abundance of Gemmatimonas and was negatively correlated with the abundance of Sphingomonas ( p < 0.05); the available nitrogen content was significantly positively correlated with the abundances of Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), was positively correlated with the abundance of Lactobacillus and Sphingomonas ( p < 0.05), and was negatively correlated with the abundances of Methylotenera , Ralstonia and Flavobacterium ( p < 0.05); the available phosphorus content was significantly positively correlated with the abundances of Lactobacillus , Gemmatimonas , Bryobacter , Candidatus and Rhizomicrobium ( p < 0.01), was positively correlated with the abundance of Bacteroides ( p < 0.05), was significantly negatively correlated with the abundances of Methylotenera , Ralstonia and Sphingobium ( p < 0.01), and was negatively correlated with the abundances of Pseudomonas and Flavobacterium ( p < 0.05); the content of organic matter was significantly positively correlated with the abundances of Lactobacillus and Gemmatimonas ( p < 0.01), was positively correlated with the abundances of Bryobacter , Candidatus and Rhizomicrobium ( p < 0.05), was significantly negatively correlated with the abundances of Methylotenera , Pseudomonas , Ralstonia and Sphingobium ( p < 0.01), and was negatively correlated with the abundance of Flavobacterium ( p < 0.05) (Table ). The results of the db-RDA showed the effect of the bacterial community on the soil bio-chemical factors in the F. taipaiensis rhizospheric soils with different cultivation years explained 68.36% of the total characteristic value (Fig. A). Further analysis showed that the pH (R 2 = 0.9005, p = 0.001), the urease activity (R 2 = 0.7792, p = 0.001), the available phosphorus content (R 2 = 0.7890, p = 0.001) and the organic matter content (R 2 = 0.8564, p = 0.001) were the dominate factors that affected by the change of the rhizospheric soil bacterial community structure, among which the pH was the most significant factor (Supplementary Table ). The average distance from the point to the center in each group of the PCoA ranking result based on the Bray–curtis distance matrix could represent the beta diversity within group. And the beta diversity of the bacterial community within group showed an upward trend after an initial decrease with the cultivation year (Fig. B). 

The Shannon and Simpson indices are indicators reflecting the species diversity influenced by both the richness and the evenness. In this study, these two indices showed a downward trend with the increase of cultivation year which represented a decline in species diversity in the rhizospheric soil of F. taipainensis . While the Chao1 and ACE indices which reflect the species richness had no significant change with the cultivation years, the decline in bacterial species diversity could be reflected in the species evenness. Previous study has shown that agricultural planting could influence the relative abundances of some specific bacterial communities in the soil which obviously alter the evenness of the bacterial communities in the soil, while the bacterial community richness would be less affected . Another study have shown that plants can shape their rhizosperic microbiome structure by the root exudates during the growth, and then directionally change the rhizospheric microbial community structure, as a result, the microbial community structures in the planting soils gradually converged . The downward trend showed by the intra group beta diversity of the bacterial communities from sample Y1 to sample Y4 could confirm this phenomenon. The results of UPGMA, NMDS and Adonis showed that the bacterial community structures in the soils with the adjacent cultivation years were closer to each other in clustering, which indicated that the bacterial community structure in F. taipaiensis rhizospheric soil changed gradually rather than suddenly with the cultivation year. Long-term single planting is considered as a process of pathogen accumulation and soil nutrient status deterioration . The dominant bacterial genera in Y1 were mainly the species that could degrade organic carbon to promote the organic cycling in soil or participate in soil nutrient cycling and the species could produce secondary metabolites to protect plants from diseases or promote the growth of plants: Haliangium produces haliangicin that can inhibit the growth of a wide spectrum of fungal pathogen ; Ferruginibacter and Chthoniobacter play important roles in soil carbon mineralization and carbohydrate metabolism , ; species in genus Piscinibacter are methylotrophic bacteria that participate in the soil carbon cycle to increase organic matter in soil ; Germmatimonas is one of the most dominant groups in agricultural ecosystems that take part in nitrogen metabolism ; Opitutus is strongly positive correlate with NO 3− -N in the soil ; Arenimonas species have catalytic activities of phosphatase and multiple lipases ; Bacteroides has the capability to produce acid that can dissolve minerals and provide nutrients for plants . While with the increase of cultivation year, the bacterial species in rhizosphere of F. taipaiensis then gradually reflected the pathogen accumulation and soil nutrient status deterioration, for Y2 and Y3: Pseudarthrobacter , Sphingomonas and Steroidobacter are significantly positively correlated to the disease incidence of continuous cropping , Acidobacteriaceae_(Subgroup_1) is negatively correlated to the content of soil organic carbon ; for Y4 and Y5: Flavobacterium , Duganella and Pseudomonas are pathogen suppressive, and are more abundant in fungal pathogen diseased plant rhizosphere , the abundance of Acidibacter is positively correlated with the damage of plant cell . Previous studies have confirmed that pH is one of the most important factors affecting the soil bacterial community structure – . Here, our study proves that pH is also the most significant factor affected by the soil bacterial community structure and its explanation degree exceeds 90%. However, previous studies have shown that, with the increase of continuous cropping years, the pH in the planting soil would gradually decrease, resulting in soil acidification – . Therefore, the increase of the pH in the rhizosphere of F. taipaiensis with the cultivation year in this study provided a special case. Among the genera with the top 15 relative abundances, Flavobacterium and Ralstonia were significantly positively correlated with pH while Bryobacter , Candidatus and Lactobacillus were negatively correlated with pH in this study. Flavobacterium is a kind of important soil saprophytic bacteria, the ammonia produced by the process of its organic matter degradation can alkalize the soil , the increasing trend of pH from Y1 to Y4 may also associated to the increasing of Flavobacterium . In addition, the ammonia produced by the process of Flavobacterium ’s organic matter degradation may be toxic to F. taipaiensis . The suitable soil environment for the growth of Fritillaria plants is weakly acidic (pH = 6.66) , therefore, in the actual planting, the rising soil pH may cause the growth and quality decline of F. taipaiensis . The metabolic process of Ralstonia can produce alkaline substances that could rise the soil pH . The relative abundance of Ralstonia in sample Y4 was obviously higher than those in Y1, Y2 and Y3, and The pH of Y4 was 7.47, the weak alkaline environment can also associated to the increasing of Ralstonia , which can badly affect the growth of F. taipaiensis . Bryobacter , Candidatus and RB41 are affiliated to Acidobacteria, the abundances of bacteria in Acidobacteria are generally negatively correlated with soil pH , , the decline on the relative abundance of these 3 genera may also lead the soil pH rising. In addition, the relative abundance of Lactobacillus which can produce lactic acid, gradually declined to 0% in Y4 and Y5, that process hastened the rising of soil pH. Some studies found that with the cultivation years’ increase, the abundances of bacteria in the cultivated soil would decline while the abundance of fungi would rise, the soil microbial population type would turn to fungal type from bacterial type, and the fungal plant pathogens would accumulate in soil during this process – . Previous studies showed that when the relative abundance of plant pathogenic fungi Fusarium rose, the the growth of Gemmatimonas could be seriously inhibited, and the relative abundance of Pseudomonas would rise . Therefore, in this study, the decline of the Gemmatimonas ’s relative abundance and the increase of the Pseudomonas ’s relative abundance may corresponded to the accumulation of pathogenic fungi with the cultivation years. In addition, Gemmatimonas was positively correlated with the activity of phosphatase and the contents of available phosphorus and organic matter in this study. Previous study have found that Gemmatimonas can dissolve insolubilized phosphate and convert it into solubilized phosphate for plant growth , its decline in relative abundance would hasten the impoverishment of planting soil. Moreover, the bacteria that declined in the relative abundance include various kind of nitrogen fixing bacteria, potash dissolving bacteria and phosphate dissolving bacteria . Thus, the above conclusion explained the reasons for the decline of the contents of available nitrogen, available phosphorus and available potassium. Urease and protease are the key enzymes involved in the transformation of soil organic nitrogen to inorganic nitrogen, which play an extremely important role in the process of soil nitrogen cycling . Urease can promote the hydrolysis of urea to ammonia in soil, its activity is often used to characterize the intensity of soil nitrogen transformation . Protease can promote the hydrolysis of proteins and peptides to amino acids, and its activity is an important characterization of soil nitrogen supply capacity . The significant higher protease activities in Y4 and Y5 compared with those in Y1, Y2, and Y3 may corresponded to the increase of Collimonas which was positively correlated with protease activities in Y3, for the Collimonas ’s metabolites had high protease activity . The changing trend of urease activity was opposite to that of protease and the urease activity was also negatively correlated with the protease activity in this study, which was consist with the result of Chang et al . Previous study found that high protease activity in the soil could significantly increase soil total nitrogen content, while the excessive total nitrogen content could inhibit the secretion of urease . Another study showed that the soil urease activity was positively correlated with the organic matter content and was negatively correlated with the invertase activity , which was consist with our study. Therefore, in our study, the urease activity showed a changing trend similar with the organic matter content and contrary to the invertase activity. Moreover, in the actual production, the cultivation of F. taipaiensis was mainly began as seed sowing, which usually took at least 5 years from germination to harvest . And because of this long-term single planting pattern, the plant would release large amount of secondary metabolites such as phenolic acids into the soil which would be toxic to the plants and affect the soil microbial structure – . Study have confirmed that phenolic acids could reduce the number of beneficial microbial population under the condition of long-term single planting . Previous study on F. taipaiensis also found that the number of culturable beneficial bacteria in the rhizospheric soil would significantly decline in the 3rd and 4th cultivation years during the single planting , this phenomenon could be explained by the accumulation of phenolic acids that also badly affect the growth of F. taipaiensis. 

In summary, long-term single planting pattern cause the obvious change of the bacterial community structure in the F. taipaiensis rhizospheric soil, the evenness of the soil bacterial community declined, the relative abundance of the beneficial bacteria decreased and the pathogenic microorganisms accumulated. The decrease of the available nitrogen, available phosphorus and organic matters impoverished the soil. This study can provide a preliminary scientific basis for the soil management in the sustainable cultivation of F. taipaiensis , solving the germplasm decline problem appeared in the production of F. taipaiensis by combining with reasonable rotation, application of organic fertilizer, development of biological fertilizer and control of harmful pathogens. In our future research, we are going to use metagenomic analysis to identify more enzymes and proteins that may be involved in the changes on the bio-chemical properties of F. taipaiensis rhizospheric soil, moreover, we will analyze the relationship between the rhizospheric microorganisms and the growth and the medicinal quality of F. taipaiensis to find the potentially beneficial microorganisms that can be used in the sustainable cultivation of F. taipaiensis . 

Supplementary Information. 

