Aerobic hydrocarbon degradation medium (AHDM), which was used for enriching crude oil–degrading bacterial populations, contains 10 g crude oil, 10 g (NH₄)₂SO₄, 1.1 g KCl, 1.1 g NaCl, 1.97 g Na₂HPO₄, 0.22 g KH₂PO₄, 0.5 g MgSO₄.7H₂O, and 0.5 mL TES per liter. Trace element solution (TES) contains 0.56 g FeSO₄, 0.17 g MnSO₄, 0.25 g CuSO₄, 0.24 g CaCl₂, and 0.29 g ZnSO₄, per liter of the solution. Petroleum was obtained from the Xinjiang oil field in China. All other chemicals and reagents used were of analytical grade.

Soil samples polluted with crude oil were collected near the Luliang oil field in Xinjiang, China. The soil samples showed that approximately 70% of the crude oil degraded after 2–3 months in a natural environment. To enrich the crude oil–degrading bacterial consortia, 5 g of soil was cultivated for 0.5 h at 200 rpm in a 100 mL flask containing 45 mL 0.9% NaCl under aerobic conditions, at 37°C, then allowed to stand for 1 h to obtain the leachate. The leachate was transferred to the shake flasks containing 100 mL AHDM and the culture was generation 1 (G1) culture. AHDM without any bacteria was used as the control. The AHDM was autoclaved at 121°C and 30 min before the cultivation experiment. The G1 cultures was incubated on an electrical shaker at 150 rpm and 37°C for 10d. At the end of cultivation, G1 culture was continuously transferred (10% v: v) into fresh 100 mL AHDM every 10 days. The subsequent culture generations were referred to by their consecutive batch culture numbers (G2, G3, and G4).

TPH content was assessed using a gravimetric method (Capelli et al., 2001). Briefly, residual oil from the culture (100 mL) was recovered by adding 40 mL carbon tetrachloride. The lower organic phase was filtered and dehydrated with anhydrous sodium sulfate. The TPH content was quantified gravimetrically, after solvent evaporation.

The residual oil was extracted and divided into saturates, aromatics, resins, and asphaltenes, based on the Chinese National Standard SY/T 5119–2016. Briefly, n-hexane was gradually added to the residual oil. The oil sample was left to stand for more than 12 h to fully precipitate the asphaltenes and measured. The filtrate was then subjected to chromatographic separation using a silica gel-neutral alumina column. Successive elution was performed using distinct organic solvents to sequentially recover the saturated hydrocarbon, aromatic hydrocarbon, and resin fractions. Saturated hydrocarbons were eluted with n-hexane, while a mixture of dichloromethane and n-hexane was employed for eluting aromatic hydrocarbons. Resin components were separated using anhydrous ethanol and chloroform as eluents. Following complete solvent evaporation, the weights of these three fractions were individually recorded.

The medium was filtered through a sterile 0.45 μm organic membrane filter (Tianjin Jinteng experimental equipment Co. Ltd., Tianjin, China). Filtered membrane was utilized for DNA extraction. DNA extraction, PCR amplification, and sequencing were performed by Guangdong Magigene Biotechnology Co., Ltd. (Guangzhou, China). The V4 and V5 hypervariable regions of microbial 16S rRNA were amplified using the primers 515F (5′–GTGCCAGCMGCCG CGGTAA–3′) and 909R (5′–CCCCGYCAATTCMTTTRAGT–3′). The constructed libraries were pooled and sequenced using 250PE (paired–end) sequencing on an Illumina novaseq platform (Illumina, USA). Analysis of bacterial consortium composition and diversity was performed in triplicate. Demultiplexing, quality filtering, clustering into Amplicon-sequence variant (ASV), and construction of the ASV table were performed in QIIME2, v.2020.2 (Bolyen et al., 2019). The analysis yielded 1,503,255 high–quality reads (94.7% of reads averaged ≥ Q30 scores) distributed across 12 samples, with the minimum and maximum number of reads per sample being 88,513 and 171,446, respectively. The deblurring algorithm was used to construct the ASV table (Amir et al., 2017). The final ASV abundance table was then rarefied at 88513 sequences per sample through the alpha-rarefaction subcommand in QIIME2. Rarefaction curves based on observed ASV and the Good’s coverage index were used to evaluate the adequacy of the sampling depth, which was established at 88513 quality-filtered reads per sample. The feature classifier script implemented in QIIME2 was employed for the taxonomic assessment using the SILVA reference database, v. 138 (Quast et al., 2013). Alpha diversity was assessed using diversity (Shannon) and evenness indices. Beta diversity was assessed by computing the Bray–Curtis distance between the samples and by Principal coordinate analysis (PCoA). Core OTUs within G1-G4 were discerned based on specific criteria: OTUs with high frequency, present in over 80% of samples, and abundant OTUs with relative abundances exceeding 0.2% across the entirety of samples (Wu et al., 2019). The neighbor-joining phylogenetic tree of the 16S rRNA gene was constructed using MEGA11 software, employing representative sequences of highly abundant bacteria (relative abundances exceeding 1%). The resulting tree was then visualized using the Interactive Tree of Life (iTOL) platform.¹ In this study, the potential functions of microbial consortia within G1-G4 were predicted by using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2 (PICRUSt2, v2.5.2) and Functional Annotation of Prokaryotic Taxa (FAPROTAX, v1.2.4).

Weighted co–occurrence analysis with an ensemble of similarity measures (Bray–Curtis, Pearson, Spearman, and MI) and regression techniques (GBLM) (Faust et al., 2012) were used to generate microbial association networks (Mac Aogain et al., 2021).

Specifically, the input file containing the relative abundance of genera in G1–G4 was calculated by five diverse measures and obtained relationships (positive and negative) between genera in G1–G4. The Gephi and Cytoscape software was used to analyze Betweenness Centrality, Closeness Centrality, and Degree of nodes. Following the result of Banerjee et al. (2018), keystone species have the characteristics of low Betweenness Centrality, high Closeness Centrality, and high Degree. To obtain the scores of the genera of the keystone species, the three indices were homogenized and used Eq. 1:

An input file, which contained the relative abundance of all genera in G1–G4 and the degradation efficiency of TPH, saturates, and aromatics, calculated the relationship between the relative abundance of genera and the degradation efficiency of TPH, saturates, and aromatics by five diversity measures. Alpha diversity in different samples was compared using the Kruskal–Wallis pairwise test (Liddicoat et al., 2020). We used Permutational multivariate analysis of variance (PERMANOVA) to investigate the effects of time and four transfers on each component of beta diversity. Network metrics such as node degree, stress centrality, and betweenness centrality were calculated and visualized using Cytoscape. Statistical analysis and visualization of bacterial consortium and alpha and beta diversities were performed using the R packages reshape2, ggplot2, dplyr, VennDiagram, igraph, and vegan. The efficiency of crude oil degradation was assessed using Eq. 2

Data are expressed as mean ± standard deviation (SD) of three replicates. The Wilcoxon-Mann–Whitney test and Kruskal–Wallis test were performed for statistical analysis of values in different groups to determine the presence of significant differences.

In this study, we enriched a crude oil-degrading bacterial consortium generation 1 (G1) from contaminated soil, achieving a TPH degradation rate of 32.29%. The crude oil degradation capacity dramatically decreased with G1 being continuously transferred. In comparison to G1, G4 consortium only achieved a crude oil degradation rate of 23.14% within the same incubation time (Figure 1A). The saturates hydrocarbons were degraded at 18.72, 22.89, 17.82, and 14.16% by bacterial consortia G1–G4 (Figure 1B), this indicates a decline in the degradation capacity of saturated hydrocarbons after G2. While the G1 consortium was capable of degrading aromatics, G2, G3, and G4 consortia exhibited negligible degradation of aromatics (Figure 1C), suggesting a loss in aromatics degradation capacity after G1. The degradation of resins and asphaltenes showed no significant difference between the bio-treatment and control groups (Supplementary Figures S1A,B). These results illustrate that TPH degradation efficiency declined during successive transfer cultures.

The rarefaction curve analysis indicated that the sequencing depth in our study was adequate to capture the diversity of bacterial communities in each sample, ensuring representative coverage (Supplementary Figure S2A). Furthermore, statistical analysis using the Kruskal-Wallis test revealed no significant differences in the Shannon diversity index and evenness index among the G1-G4 consortia (p > 0.1; Supplementary Figures S2B,C). When assessing beta diversity based on the Bray-Curtis distance, PERMANOVA analysis demonstrated no significant difference (p > 0.1) in bacterial community structure between G1-G4 consortia (Supplementary Figure S2D).

Fourteen shared OTUs were observed in the four generations, while 24, 10, 9, and 16 unique OTUs were found, respectively, in G1, G2, G3, and G4 (Figure 2A). There were 5 OTUs in 14 shared OTUs identified as the regional core OTUs, accounting for 20.0 to 51.6% of the bacterial abundance of four generations (Figure 2B). All the 59 unique OTUs were classified to 5 phyla, including Actinobacteriota, Bacteroidota, Chloroflexi, Firmicutes, and Proteobacteria (Figure 2C). Compared to another generations, the unique OTUs of G1 were predominantly concentrated in Firmicutes and Proteobacteria. Those of G2 were mainly centered around Bacteroidota, and G4 primarily exhibited a concentration of Actinobacteriota.

In Figure 2D, the representative sequences of highly abundant bacteria (relative abundances exceeding 1%) in this study were selected to construct the phylogenetic tree. We found that the bacterial consortia of G1-G4 had the characteristic of high diversity. In G1, dominant OTUs were primarily concentrated in the Bacillus, Dietzia, Novosphingobium, Xanthomonas genus, and the Pseudomonadales order with Acinetobacter and Pseudomonas genera. In G2, predominant OTUs were mainly centered around Pseudomonas, Dietzia, Novosphingobium, and Paenibacillus. For G3, dominant OTUs were primarily found in Novosphingobium, Sporosarcina genus, and the Enterobacterales order. In G4, predominant OTUs were focused on the Luteibacter genus and the Enterobacterales order.

The analysis of the bacterial community revealed that Proteobacteria (59.35% ± 17.64%), Firmicutes (25.29% ± 9.18%), and Actinobacteriota (14.25% ± 8.63%) were the predominant phyla in G1-G4 consortia (Figure 2E). The dominant genera exhibited variations among the consortia. In G1, prevalent genera included Pseudomonas (18%), Bacillus (27.9%), Dietzia (20%), and Novosphingobium (8.27%). However, in G4, the dominant genera shifted to Serratia (22%), Luteibacter (18.9%), Sporosarcina (6.97%), Bacillus (4.77%), Pseudomonas (7.93%), and Novosphingobium (6.89%) (Figure 2F). Statistical analysis using the Wilcoxon-Mann–Whitney test indicated significant differences in the relative abundance of certain genera between G4 and the earlier generations. Specifically, the relative abundance of Pseudomonas in G4 was significantly lower than in G1 and G2, while Bacillus in G3 and G4 showed significantly lower abundance compared to G1 and G2. Dietzia in G4 exhibited a significantly lower abundance compared to G1-G3, and Novosphingobium in G4 showed a significantly lower abundance than in G3. By contrast, Serratia and Luteibacter in G4 had significantly higher abundance than in G3, and Sporosarcina in G4 showed significantly higher abundance than in G1 (Figure 2F).

FAPROTAX and PICRUSt2 software were be used to predict the function of bacterial communities. FAPROTAX was used to predict the biochemical cycle of environmental samples. To display potential function of crude-oil degradation of bacterial communities, C transformation functions were promoted, including chemoheterotrophic, aerobic chemoheterotrophic, and aromatic_compound degradation processes (Figure 3). Our findings demonstrated that as successive generations progress, the aerobic chemoheterotrophic process of G1 and G3 became significantly more active in G4. The aromatic degradation processes of G1 exhibited pronounced activation in G2, G3, and G4. Notably, no significant differences in chemoorganotrophic were observed across G1 to G4.

Alkanes primarily undergo intricate transformations facilitated by specific enzymatic systems such as alkane hydroxylases, ethanol dehydrogenases, acetaldehyde dehydrogenases, and acetyl-CoA. These processes culminate in their eventual entry into the fatty acid beta-oxidation pathway, leading to their complete oxidation into carbon dioxide and water (Rojo, 2009). The degradation of polycyclic aromatic hydrocarbons (PAHs) involves a series of sequential reactions including hydroxylation, dehydrogenation, isomerization, and ring cleavage. This intricate sequence results in the generation of intermediates that become integrated into the tricarboxylic acid (TCA) cycle. Subsequently, these intermediates undergo enzymatic conversions orchestrated by microorganisms, ultimately leading to the production of carbon dioxide and water. Notably, the ring cleavage process constitutes a pivotal step in the degradation of PAHs, involving two distinct pathways (Varjani, 2017): (1) ortho-cleavage catalyzed by catechol 1,2-dioxygenase (C12O) at ortho positions, and (2) meta-cleavage catalyzed by catechol 2,3-dioxygenase (C23O) at meta positions. Utilizing the PICRUSt2 software, we acquired the abundance of genes encoding these enzymes. As illustrated in the Figure 4, the relative abundance of the alkane hydroxylase gene (K00496) in G1 notably surpassed that in G3 and G4, while both G2 and G3 exceed G4. The acetaldehyde dehydrogenase gene (K00001) in G1 and G3 also exhibited significantly higher abundance than in G4. Furthermore, the acetyl-CoA gene (K01897), along with the C12O gene (K03381) and C23O gene (K07104) involved in aromatic hydrocarbon metabolism, exhibits significantly higher abundance in G1, G2, and G3 compared to G4. Notably, the ethanol dehydrogenase gene (K00128) displays negligible variation between G1 and G4. Consequently, the trend observed across G1 to G4 generations unveils a marked decline in the relative prevalence of genes linked to alkane and aromatic hydrocarbon metabolism, excluding the ethanol dehydrogenase gene.

Based on key-score values and the relationship between genera, we constructed Figure 5. Pseudomonas exhibited the highest key-score in G1, followed by Arthrobacter in G2, Taonella in G3, and Lysobacter in G4 (Figures 5A–D). Their relative abundances were 18, 0.21, 0.03, and 1.77%, respectively (Supplementary Table S1). These findings indicate that the presence of key functional species in the consortium is not solely dependent on their relative abundance. Alternatively, important functional species can exist as either dominant populations with high relative abundance or non-dominant populations with low relative abundance.

Table 1 showed the correlation between the genus of top key scores in G1-G4 and the degradation efficiency of TPH, saturates, and aromatics. The correlations ranged from 15.42 to −9.54 for TPH, 11.91 to −9.79 for saturates, and 8.03 to −5.96 for aromatics. These results suggest a decline in the degradation capacity of TPH, saturated hydrocarbons, and aromatic hydrocarbons by the dominant species in the consortia. In Table 2, the highest positive correlation with TPH and saturated hydrocarbons degradation was observed for the Dietzia genus, with correlation values of 19.48 and 18.38, respectively (p < 0.05). Similarly, the Bacillus genus showed the highest positive correlation with aromatic hydrocarbons degradation, with a correlation value of 16.56 (p < 0.05). Conversely, in Table 3, Serratia exhibited the highest negative correlation with TPH and saturated hydrocarbon degradation, with correlation values of −19.38 and − 15.78, respectively (p < 0.05). Additionally, Psychrobacter showed the highest negative correlation with aromatic hydrocarbon degradation, with a correlation value of −8.92 (p < 0.05). These results emphasize the important roles of the Dietzia and Bacillus genera in the degradation of TPH, saturated hydrocarbons, and aromatic hydrocarbons, while highlighting the negative influence of Serratia and Psychrobacter on the degradation process.

The key scores of the Dietzia genus decreased from 3.62E-02 to 2.58E-02 during successive transfer cultures (Table 4). In contrast, the key scores of the Serratia and Psychrobacter genera increased from 0 to 2.87E-02 and 3.72E-02, respectively, after G2 (Table 4). The decline in the degradation of TPH and saturated hydrocarbons can be ascribed to the gradual weakening of the relationship between Dietzia and other genera in the consortia, as well as the progressive strengthening of the relationship between Serratia and Psychrobacter with other genera. These changes in inter-genera relationships within the consortia have significant implications for the overall degradation efficiency of TPH and saturated hydrocarbons.

The Bacillus genus, exhibiting the highest positive correlation (16.56) with aromatic hydrocarbon degradation, exhibit positive correlation with Mycetocola genus only in G1 (Figures 5A–D). This correlation pattern aligns with the observed degradation pattern of aromatic hydrocarbons during successive transfer cultures, suggesting that the relationship between these two genera plays a role in influencing the degradation of aromatic hydrocarbons.