Based on the relative abundance table of different classification levels, the composition and relative abundance of gut bacteria in 6-OHDA-induced PD rat (6-OHDA), 6-OHDA-induced PD rat treated with 5 mg/kg levodopa (6-OHDA+Ldopa5), 6-OHDA-induced PD rat treated with 10 mg/kg levodopa (6-OHDA+Ldopa10) and control groups were examined focusing on the 10 most abundant taxa at phylum-, class-, order-, family-, genus-and species-levels, respectively, and the rest taxa were set as others. The relative abundance histogram of the corresponding taxa annotation results of each group was drawn at different classification levels (Figure 1). In detail, at the phylum level, the relative abundance of Fibrobacteres increased apparently in 6-OHDA+Ldopa10 group compared to 6-OHDA group (Figure 2A). At the class level, the relative abundance of Bacilli decreased obviously in 6-OHDA+Ldopa10 group compared to 6-OHDA+Ldopa5 group (Figure 2B), and the relative abundance of Negativicutes increased strikingly in 6-OHDA+Ldopa10 group compared to 6-OHDA group (Figure 2C). At the order level, the relative abundance of Lactobacillales decreased notably in 6-OHDA+Ldopa10 group compared to 6-OHDA+Ldopa5 group (Figure 2D). At the family level, the relative abundance of Lactobacillaceae decreased prominently in 6-OHDA+Ldopa10 group compared to 6-OHDA+Ldopa5 group (Figure 2E). At the genus level, the relative abundance of Lactobacillus decreased remarkably in 6-OHDA+Ldopa10 group compared to 6-OHDA+Ldopa5 group (Figure 2F), and the relative abundance of Roseburia increased prominently in 6-OHDA+Ldopa10 group compared to the 6-OHDA+Ldopa5 group (Figure 2G). At the species level, the relative abundance of Lactobacillus johnsonii decreased significantly in 6-OHDA+Ldopa10 group compared to 6-OHDA and 6-OHDA+Ldopa5 groups (Figure 2H), while the relative abundance of Lactobacillus murinus and Prevotella sp. CAG:485 decreased distinctly in 6-OHDA+Ldopa10 group compared to 6-OHDA group (Figures 2I,J). These results indicated that an increase in dosage of levodopa could potentially change the gut bacteria composition in 6-OHDA-induced PD rats.

Next, LEfSe was performed to analyze the differentially abundant of gut bacteria among the four groups. In order to screen the species biomarkers with significant differences between groups, the rank sum test was used to detect the different species in different groups, and linear discriminant analysis (LDA) was used to achieve dimension reduction and assess the impact of the different species, namely, LDA score. LEfSe analysis results of different species between groups included histogram of LDA value distribution and evolutionary clade diagram (phylogenetic distribution), which showed that 75 bacterial communities were different among control, 6-OHDA, 6-OHDA+Ldopa5 and 6-OHDA+Ldopa10 groups. There were 30 biomarkers in control, 13 biomarkers in 6-OHDA, 10 biomarkers in 6-OHDA+Ldopa5 and 19 biomarkers in 6-OHDA+Ldopa10 (Figure 3A). The biomarkers in control group mainly included Bacteroides_sp_CAG_927 (species), Clostridium (genus), Lactobacillus_sp_ASF360 (species), Unclassified_Bacteroidales (family), and Romboutsia_ilealis (species), among which the most abundant was Bacteroides_sp_CAG_927 (species). The biomarkers in 6-OHDA group mainly included Lactobacillus_murinus (species), Firmicutes_bacterium_CAG_424 (species), Faecalibacterium (genus), Firmicutes_bacterium_CAG_646 (species), and Alistipes (genus), among which the most abundant was Lactobacillus_murinus (species). The biomarkers in 6-OHDA+Ldopa10 group mainly included Prevotella (genus), Prevotellaceae (family), Roseburia_sp_CAG_309 (species), Prevotella _sp_P2_180 (species), and Prevotella _sp_P5_92 (species), in which the most abundant was Prevotella (genus). The biomarkers in 6-OHDA+Ldopa5 group mainly included Lactobacillus (genus), Lactobacillus_johnsonii (species), Bacteroides_sartorii (species), Lactobacillus_animalis (species), and Bifidobacteriates (order), in which the most abundant was Lactobacillus (genus). The cladogram showed significant changes in each group (Figure 3B).

Resistance genes are ubiquitous in both human gut microbes and other environmental microbes. Comprehensive Antibiotic Resistance Database (CARD) is a newly emerged resistance gene database in recent years. The core component of CARD is ARO, which integrates the antibiotic sequence, antibiotic resistance and function mechanism, et al. Using the Resistance Gene Identifier program, the unigenes were aligned to the CARD. Based on the relative abundance of resistance gene, the common top 20 AROs identified in 6-OHDA, 6-OHDA+Ldopa5, 6-OHDA+Ldopa10 and control groups were tetQ, tetW/N/W, tetO, adeF, AAC6-Ib-Suzhou, OXA-397, tetW, tet44, tet32, vanTG, sul3, ErmB, tet40, Escherichia_coli_acrR, adeR, mdtB, MCR-5, tetBP, ErmF and MexD (Figure 4A). In comparison to the 6-OHDA group, the relative percentage of AAC6-Ib-Suzhou in the 6-OHDA+Ldopa10 group was greater (Figure 4B). In comparison to the 6-OHDA+Ldopa5 group, the relative percentage of tetW in the 6-OHDA+Ldopa10 group was lower (Figure 4C). In comparison to the 6-OHDA group, the relative percentage of vanTG in the 6-OHDA+Ldopa10 group was lower (Figure 4D). To more intuitively observe the proportion of ARO abundance in each group on the whole and more intuitively display the overall distribution of ARO abundance, the AROs with the top10 abundance was selected to draw the overview circle diagram (Figure 4E).

In order to investigate the differences in the number of ARG and ARO between the sample groups, the box maps of the number of ARG and ARO between the sample groups were drawn, which showed that the number of ARG and ARO among all groups did not significantly differ (Figures 5A,B). A Venn diagram was also used to examine how the four groups’ shared ARGs were composed. The Venn diagram for the ARG number classified by the four component approaches is depicted in Figure 5C. The four component methods complemented each other in classification. Four, 7, 14 and 6 ARGs were uniquely correctly classified by control, 6-OHDA, 6-OHDA+Ldopa5 and 6-OHDA+Ldopa10 groups, respectively. Moreover, the Venn diagram for the ARG number classified by the three component approaches is depicted in Figure 5D. The three component methods complemented each other in classification. Eleven, 13 and 26 ARGs were uniquely correctly classified by control, 6-OHDA and 6-OHDA+Ldopa5 groups, respectively. In addition, the Venn diagram for the ARG number classified by the three component methods. The three component methods complemented each other in classification. Twelve, 14 and 18 ARGs were uniquely correctly classified by control, 6-OHDA and 6-OHDA+Ldopa10 groups, respectively (Figure 5E). Furtherly, the Venn diagram for the ARG number classified by the two component methods. The two component methods complemented each other in classification. Thirty-three and 26 ARGs were uniquely correctly classified by control and 6-OHDA groups (Figure 5F).

We housed 6-week-old male SD rats with free access to food and water that weighed 180–200 g. The institutional Animal Care and Use Committee of Zunyi Medical University (Zunyi, China) gave its approval to all animal experiments.

The experiments were carried out as formerly reported (Zheng et al., 2021). In order to target the unilateral SNpc (coordinates AP-5.2 mm, ML-2.1 mm, DV-8.0 mm) relative to the bregma, rats were anesthetized and then immobilized in a stereotaxic frame, 6-OHDA solution was injected over 3 min with an infusion rate of 1 μl/min, followed by 3 min of equilibrium before retracting the needle. Prior to the experiment, all rats had recovered after 3 weeks. Then, for 12 weeks, a daily intraperitoneal injection of 5 mg/kg or 10 mg/kg levodopa was administered.

Following a 12-week course of levodopa therapy, feces were collected according to the previously reported methodology (Zhang and Wang, 2017), Each sample was briefly lysed using lysozyme-containing CTAB buffer. Chloroform, isoamyl alcohol (24:1), and isopropanol were used to extract the whole DNA. The precipitate was then twice rinsed with 75% ethanol before being dissolved in ddH₂O, RNase A was then used to degrade the RNA.

For sample testing, on 1% agarose gels, the level of DNA degradation and potential contamination were observed, DNA concentration was measured using Qubit® dsDNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, United States), DNA concentrations greater than 1 ug are used to construct libraries with an OD value of 1.8 to 2.0. The input material for the DNA sample preparations was a total of 1 g of DNA per sample. Utilizing the NEBNext® Ultra™ DNA Library Prep Kit for Illumina (NEB, United States) in accordance with the manufacturer’s instructions, sequencing libraries were created, and index codes were added to assign sequences to specific samples. The DNA sample was broken up to 350 bp-sized fragments using sonication, after which it was end-polished, A-tailed, and ligated with the full-length adaptor for Illumina sequencing and additional PCR amplification. Finally, PCR products were purified using the AMPure XP system, and libraries underwent size distribution analysis with an Agilent 2,100 Bioanalyzer and real-time PCR quantification. On a cBot Cluster Generation System, the index-coded samples were clustered in accordance with the manufacturer’s recommendations. The library preparations were sequenced on an Illumina HiSeq platform after cluster generation, and paired-end reads were produced.

To obtain the clean data needed for further analysis, Readfq (V8¹) was used to preprocess the raw data collected from the Illumina HiSeq sequencing technology. These were the precise processing steps: reads were first removed if they contained low-quality bases (default quality threshold value ≤38) above a predetermined portion (default length of 40 bp); then, reads were removed if the N base reached a predetermined percentage (default length of 10 bp); and finally, reads were removed if they shared the overlap with the adapter above a predetermined length (default length of 15 bp). Given the potential for host contamination in samples, clean data had to be sent to the host database, and Bowtie2.2.4 software was used by default to filter out reads of host origin (Bowtie2.2.4²), the following were the parameters: –end-to-end, −-sensitive, -I 200, −X 400 (Karlsson et al., 2012, 2013).

The SOAPdenovo software (V2.04³) was used to assemble and analyze the clean data(Luo et al., 2012), the following were the parameters: -d 1, −M 3, −R, −u, -F, -K 55 (Scher et al., 2013; Qin et al., 2014; Brum et al., 2015; Feng et al., 2015). MEGAHIT software (v1.0.4-beta) could be used to assemble the Clean Data for samples taken from complex environments, such as water, soil, and so forth, and the parameters are-presets meta-large (−- min-count 2 --k-min 27 --k-max 87 --k-step 10). Then, the built Scaftigs were severed from N connection and the Scaftigs were left without N (Mende et al., 2012; Nielsen et al., 2014; Qin et al., 2014). All samples’ clean data were compared to each Scaffolds correspondingly using Bowtie2.2.4 software to acquire the PE reads not used and the parameters were: --end-to-end, −-sensitive, -I 200, −X 400 (Qin et al., 2014).

The MetaGeneMark (V2.10) software, available at http://topaz.gatech.edu/GeneMark/, predicted all ORFs for the scaftigs (≥500 bp) assembled from single and mixed samples, and filtered with the length information shorter than 100 nt from the anticipated result via default parameters (Qin et al., 2010, 2014; Li et al., 2014; Nielsen et al., 2014; Zeller et al., 2014). For ORF predicted, CD-HIT software (V4.5.8⁴; Li and Godzik, 2006; Fu et al., 2012) was adopted to obtain the unique initial gene catalog [the genes here referred to the nucleotide sequences coded by unique and continuous genes (Sunagawa et al., 2015)]. The parameters option (Zeller et al., 2014; Sunagawa et al., 2015) were-c 0.95, −G 0, -aS 0.9, −g 1, −d 0. Using Bowtie2.2.4, the clean data from each sample was mapped to the initial gene catalog. The number of reads was changed to reflect the number of genes mapped in each sample using the parameter settings --end-to-end, −-sensitive, -I 200, −X 400 (Li et al., 2014; Qin et al., 2014). The gene which the number of reads ≤2 (Qin et al., 2012; Li et al., 2014) was filtered in each sample and the gene catalog (Unigenes) was eventually used for subsequent analysis. Based on the amount of mapped reads and the length of gene, statistic analysis on the abundance information of each gene in each sample was done. The format was as follow, “r” represented the number of reads mapped to the genes and “L” represented gene’s length (Qin et al., 2010; Karlsson et al., 2012; Cotillard et al., 2013; Le Chatelier et al., 2013; Zeller et al., 2014; Villar et al., 2015). The abundance of each gene in each sample in the gene catalog served as the basis for the fundamental data of statistics, core-pan gene analysis, correlation analysis of samples, and Venn diagram analysis of number of genes.

The DIAMOND program (V0.9.9⁵; Buchfink et al., 2015) was used to blast the Unigenes to the sequences of Bacteria, Fungi, Archaea, and Viruses from the NCBI NR database (Version: 2018-01-02⁶) using the parameter-e 1e-5. The LCA algorithm was used for system classification of MEGAN (Huson et al., 2011) software to ensure the species annotation information of sequences. For the finally aligned results of each sequence, each sequence may have multiple aligned results, and the result of which the e value the smallest e value * 10 (Oh et al., 2014) was chosen. Based on the LCA annotation result and the gene abundance table, a table including the number of genes and the abundance information of each sample in each taxonomy hierarchy (kingdom, phylum, class, order, family, genus, and species) was created. The abundance of one specie in one sample was labeled as the sum of the gene abundances for the specie. A specie’s gene number in a sample was equal to the number of genes with nonzero abundance. Based on the abundance table of each taxonomic hierarchy, Krona analysis was used to display the generation situation of relative abundance, the abundance cluster heat map, PCA (Avershina et al., 2013; R ade4 package, Version 2.15.3), and NMDS (Noval Rivas et al., 2013; R vegan package, Version 2.15.3). Anosim analysis was used to examine the difference between groups (R vegan package, Version 2.15.3). To find the distinct species among groups, metastats and LEfSe analysis were utilized. In Metastats analysis, the permutation test across groups was performed to obtain p value for each taxonomy, and then the Benjamini and Hochberg False Discovery Rate was calculated to correct p value and obtain q value (White et al., 2009). The LEfSe program carried out the LEfSe analysis (the default LDA score is 3) (Segata et al., 2011). Finally, a random forest of 6-OHDA was created using the random forest (RandoForest) (Breiman, 2001; R pROC and randomForest packages, Version 2.15.3). By using MeanDecreaseAccuracy and MeanDecreaseGin, the relevant species were eliminated. Each group was then cross-validated (by default, 10 times), and the ROC curve was generated.

To blast Unigenes to a functioning database, DIAMOND software (V0.9.9) was used using the parameter setting of BLASTp, −e 1e-5 (Li et al., 2014; Feng et al., 2015). Functional databases were exempted from the study, including the KEGG (Kanehisa et al., 2006, 2014) database (Version 2018-01-01⁷), the eggnog (Powell et al., 2014) database (Version 4.5⁸), and the CAZy (Cantarel et al., 2009) database (Version 201,801⁹). The best Blast Hit for each sequence’s blast result was utilized for further investigation (Li et al., 2014; Bäckhed et al., 2015; Feng et al., 2015). A statistical study of the relative abundance of several functional hierarchies was annotated, with the relative abundance of each functional hierarchy being equal to the total of the relative abundances. The gene number table for each sample in each taxonomic hierarchy was obtained based on the function annotation result and gene abundance table. A function’s gene number in a sample matched the gene number that was ascribed to it, and the abundance was non-zero. The number of annotated genes was counted, the overall relative abundance situation was displayed, an abundance cluster heat map was displayed, and decrease-dimension PCA and NMDS analyzes were performed using the abundance tables of each taxonomic hierarchy. Anosim study of the difference between groups (inside) based on functional abundance, comparative analysis of metabolic pathways, and analyzes of functional difference between groups using Metastats and LEfSe were also carried out.

With the parameter setting of BLASTp, evalue ≤1e-30, the Unigenes database¹⁰ (McArthur et al., 2013; Martínez et al., 2015; Jia et al., 2017) was aligned to the CARD database using the Resistance Gene Identifier (RGI) program. The aligned result was used to count the relative abundance of ARO. The abundance bar charts, abundance cluster heatmap, and the numerical difference between groups of resistance genes were all shown using the abundance of ARO. The distribution of resistance genes’ abundance in each sample, an analysis of the resistance genes’ species attribution, and a study of the resistance mechanisms of the resistance genes were also carried out.

The data were shown as mean ± SEM. For the statistical comparison, one-way analysis of variance (ANOVA) was performed using SPSS statistical software. The Dunnett’s test was applied to all pairwise comparisons of means where analysis of variance revealed significant differences. Statistical significance was defined as p < 0.05.