Genomic alterations in Bacteroides fragilis favor adaptation in colorectal cancer microenvironment

Abstract

Background

The occurrence and development of colorectal cancer (CRC) is an incredibly long process that involves continuous changes in the tumor microenvironment. These constant changes may ultimately result in genetic alterations and changes in the metabolic processes of some symbiotic bacteria as a way to adapt to the changing environment. Patients with CRC exhibit an altered abundance of Bacteroides fragilis (B. fragilis) as indicated by several studies. To better understand the genomic characteristics and virulence spectrum of B. fragilis strains in tumor tissues, B. fragilis strains were isolated from tumor and paracancerous tissues of CRC patients.

Methods

The isolates were identified using 16 S rRNA sequencing, morphological analysis, physiological and biochemical characterization and PCR, and they were then subjected to whole genome sequencing (WGS) analysis.

Results

A strain of B. fragilis enterotoxin (BFT) bft1-producing ZY0302 and a non-enterotoxin-producing B. fragilis ZY0804 were isolated from cancerous and paraneoplastic tissues, respectively. Analysis based on the core and nonessential genes showed that the genomic profiles of the isolates, ZY0302 and ZY0804, differed from those of B. fragilis from other tissue sources. This core and the co-evolution of non-essential genes may be the result of their adaptation to fluctuations in the tumor microenvironment and enhancing their survival. In addition, the ZY0302 and ZY0804 genomes underwent extensive horizontal gene transfer and varying degrees of genomic rearrangements, inversions, insertions, and deletion events, which may favor the enhancement of bacteria’s ability to adapt to environmental changes. For instance, the virulence factors, such as the capsular biosynthesis gene clusters and components of the type IV secretion system, acquired through horizontal gene transfer, may facilitated B. fragilis in evading immune responses and managing oxidative stress. Moreover, our analysis revealed that multiple virulence factors identified in the isolates were mainly involved in bacterial adhesion and colonization, oxidative stress, iron acquisition, and immune evasion. This observation is worth noting given that enzymes such as neuraminidase, lipase, hemolysin, protease, and phosphatase, along with genes responsible for LPS biosynthesis, which are recognized for their association with the virulence of B. fragilis, were prevalent among the isolates.

Conclusions

In summary, it is our assertion that the alterations observed in both core and nonessential genes of B. fragilis, which have been isolated from tissues of colorectal cancer patients, along with significant instances of horizontal gene transfer to the genome, are likely intended to enhance adaptation to the evolving conditions of the tumor microenvironment. This study may provide new insights into the interaction between B. fragilis and the CRC microenvironment.

Introduction

As one of the most common cancer types, colorectal cancer (CRC) is recognized as the third most diagnosed cancer and also ranks as the second leading cause of cancer-related deaths, highlighting its significant threat to human health [1]. Currently, the pathogenesis of CRC remains obscure as its occurrence is complex and involves multiple factors such as genetics, environment, diet and lifestyle. As the key metabolic and immunomodulatory factors of the host, intestinal microorganisms participate in crucial regulatory roles where they mediate environmental factors and the occurrence and development of CRC [2]. The development of sequencing technology has provided evidence to show that the microbial diversity in tumor tissues of patients with CRC is altered and that this change may be an important factor in promoting the development of CRC [3, 4].

B. fragilis is a commensal Gram-negative facultative anaerobic bacterium found mainly in the colon, accounting for approximately 1% of the intestinal microbiota [5]. B. fragilis is also an opportunistic pathogen that has been strongly associated with the development of acute diarrhea, inflammatory bowel disease (IBD), and CRC, and it is one of the major causative agents of anaerobic sepsis [6]. B. fragilis is classified into enterotoxin-producing B. fragilis (ETBF) and non-enterotoxin-producing B. fragilis (NTBF), depending on its ability to secrete B. fragilis enterotoxin (BFT) [5]. BFT is a ~ 20-kDa zinc metalloprotease [7], which cleaves E-calmodulin, disrupts inter-epithelial cellular barriers [8], and induces c-Myc expression and cellular proliferation [9], as well as activating helper T cell 17 (Th17) response and NF-κB signaling pathway [10], thereby promoting intestinal inflammation and colorectal carcinogenesis. In addition, ETBF induces cellular DNA damage, thereby increasing the risk of cancer development [11]. Recent studies have also shown that ETBF and its generated BFT are associated with high levels of CpG island methylation in CRC and may promote colorectal tumorigenesis through the serrated tumor pathway [12, 13]. NTBF, however, has been recognized as a potential probiotic that promotes the production of short-chain fatty acids [14] and inhibits the growth of other pathogenic microorganisms [15]. Additionally, it contributes to the diversity of intestinal flora and prevents intestinal epithelial damage, which in turn inhibits intestinal inflammation and inflammation-associated colorectal cancer [16, 17]. Polysaccharide A (PSA) is a major effector of NTBF, and TLR2 signaling plays a crucial role in mediating the protective effect against colitis-associated CRC [17].

Related studies have shown that NTBF strains isolated from different tissues may exhibit distinct biological characteristics. For example, the use of proteomics indicates that NTBF strains from patients with polyps are different from NTBF strains isolated from patients without polyps. NTBF strains from patients with polyps are rich in LPS biosynthetic genes, which may increase their ability to activate the immune system and cause inflammation [18]. Various changes may occur during the development and progression of CRC such as the growth of the transformed tissues and the secretion of various tissue factors and cellular metabolites. Consequently, the local environment may change, making it difficult for some commensal bacteria to survive. As a means of adaptation to the altered microenvironment, commensal bacteria may switch on new genetic and metabolic programs [18]. A research study found that B. fragilis possesses numerous DNA inversion systems distributed across its genome. This genetic feature allows the organism to selectively express a diverse array of cell surface structural proteins, which in turn enhances its ability to adapt to varying nutritional environments, assists in evading the host immune response, and boosts its potential for colonization and infection [19, 20]. In conclusion, B. fragilis can adapt to tissue-specific and microenvironment changes through genomic changes.

Researchers have long speculated about the role of gut bacteria in the development of CRC, and the isolation and culture methods of target bacteria have made the exploration of the etiology of CRC easier and achievable [21]. Koch’s postulates have long proposed that microbial culture is the gold standard for detecting bacteria, but most bacteria in the colon are “unculturable”, which makes it a great challenge to further explore the pathogenic mechanism of intestinal pathogens [21]. The isolation of bacteria from clinical specimens is usually affected by many factors, as such studies have pointed out that the ideal specimen (tissue or fluid) for the isolation of bacteria should be obtained during surgery [22]. When isolating and culturing B. fragilis, proper anaerobic collection containers and methods must be ensured. B. fragilis is an obligate anaerobe and requires specialized transport conditions, such as dedicated anaerobic transport media [22]. The anaerobic nature of B. fragilis requires specialized environments for its isolation and culture, which makes research challenging. Therefore, the primary goal of this study is to successfully and effectively isolate B. fragilis from CRC tumors and paracancerous tissues and culture it.

At present, B. fragilis used in CRC research is obtained from skin tissue after a traumatic infection, fecal samples, patients with septicemia, or from other species like piglets; however, isolates from tumor sites of patients are rarely reported. To narrow the gap in infrequent reports pertaining to isolates from tumor sites, the aim of this study was to isolate B. fragilis from tumors and paracancerous tissues of CRC patients. The whole genome sequences of these isolates were analyzed to investigate their evolutionary features, as well as the virulence genes of the isolates. The aim is to provide a better understanding of the genomic characteristics and virulence profiles in CRC tumor tissues.

Materials and methods

Specimen collection

Samples were collected from tumor tissues and paracancerous tissues of CRC patients from the Affiliate Hospital of Zunyi Medical University (Guizhou, China), and informed and signed consent was obtained from all patients. After collection into sterile sampling bags, the samples were stored in ice boxes and immediately sent to the laboratory for testing. This study was approved by the Medical Ethics Review Committee (ZunYi ethical examination No. [2023] 1–014).

Isolation and purification of B. fragilis

According to the method reported in literature [23, 24], the specific scheme of bacterial isolation is as follows: Tissue samples were obtained from the center of the tumor and paracancerous tissue using an aseptic surgical knife cauterized by an alcohol lamp. The tissue sections were inoculated into Brain Heart Infusion Broth (BHI Broth) medium with an aseptic inoculation loop, and the anaerobic culture was incubated at 37℃ for 18–24 h. Then the bacterial solution was inoculated into Bacteroides Bile Esculin (BBE) blood plates and anaerobically incubated at 37℃ for 24–48 h. The colony morphology, color, size and surface properties were studied and observed. At least three single colonies (slightly convex morphology, smooth and shiny surface, neat edge, brown or off-white, and black surroundings) were selected from each plate and inoculated into the CDC anaerobic blood plate. The same colony was inoculated into three separate plates and incubated under different conditions, including anaerobic, microaerobic, and aerobic environments at 37 °C for 24–48 h. Observed the growth status of the colonies in the three environments, selected the colonies that grow only in anaerobic environment, inoculated them on BBE blood plates and incubated them anaerobically for 24–48 h. After that, picked a single colony and inoculated it into 10 mL of BHI Broth medium, and then incubated it anaerobically at 37℃ for 18–24 h to get the purified strains. Blood plates used for isolation contained 5% sheep blood, vitamin K1 (0.5 mg/L), and hemin (5 mg/L). Each purified strain was made into a glycerol stocks using 0.25 mL sterile glycerol, 0.75 mL sterile PBS and 0.5 mL bacterial solution mixed thoroughly and preserved at -20℃ and liquid nitrogen.

Identification of target strains

Morphological, physiological and biochemical characteristics

After anaerobic incubation at 37℃ for 24 h in BHI Broth medium, the isolated strains were then inoculated on CDC anaerobic blood plate containing 5% sheep blood, vitamin K1 (0.5 mg/L), and hemin (5 mg/L), as well as BBE blood plates. These plates were then incubated at 37℃ for 24–48 h, and their morphological characteristics were studied and observed. The isolates were stained with Gram staining solution (HuanKai Microbial, Guangdong, China), and the shape and size of the bacteria were observed using an optical microscope (OLYMPUS BX-53, Tokyo, Japan). The physiological and biochemical characteristics of the isolated strains were evaluated using microbial microbiochemical identification tubes (HuanKai Microbial, Guangdong, China).

DNA extraction, PCR amplification and phylogenetic analysis

The isolated strains were inoculated in BHI Broth medium and anaerobically incubated at 37℃ for 18–24 h. Genomic DNA was extracted utilizing Ezup column bacteriological DNA extraction kit (Sangon Biotech, Shanghai, China). To identify the species of isolated strains, the sequences of 16 S ribosomal RNA gene (16 S rRNA) were amplified using universal primers, 27 F (5′-AGAGTTTGATCATGGCTCAG-3′) and 1492R (5’-GGTTACCTTGTTACGACTT-3’). The final PCR master mix was 50 µL, including 5 µL 10 × PCR Buffer (Mg2 + Plus), 0.5 µL r Taq (250U, 5U/mL), 2 µL dNTP Mixture (2.5 mM each), 0.5 µL each of upstream and downstream primers (10 µM), 2.5 µL genomic DNA, and 39 µL ddH2O. PCR amplification was carried out using BIO-RAD T100TM Thermal Cycler, the samples were pre-denatured at 95℃ for 5 min, denatured at 95℃ for 30 s, annealed at 55℃, extended for 1 min at 72℃ for 30 cycles, and finally extended for 10 min at 72℃. After the completion of the process, 5 µL of the PCR product was used for 1% agarose gel electrophoresis to confirm the fragment amplified by PCR. Subsequently, the PCR product was sent to Sangon Biotech (Shanghai) Co., Ltd for sequencing. The sequencing results were analyzed using BLAST alignment by the National Center for Biotechnology Information (NCBI). According to the BLAST results, 16 S rRNA sequences which were closely related to the sequences in this study were retrieved from GenBank for phylogenetic analysis. The phylogenetic tree was constructed by the maximum likelihood method by MEGA11 software, and 1000 Bootstrap repeated samples were used to evaluate the reliability and statistical significance of the constructed phylogenetic tree.

After the 16 S rRNA results confirmed that the strain was B. fragilis, the BFT, bft-1, bft-2, and bft-3 gene sequences of the strains were amplified by PCR using specific primers, and the type of the isolated strain was determined by 1% agarose gel electrophoresis. The primer sequence is shown in Table 1. The final PCR master mix and conditions employed were similar to those used to amplify the 16 S rRNA sequence.

Table 1 Primer sequences of B. fragilis toxin genotypes

Genome sequencing and assembly

Genomic DNA was extracted using the Ezup Column Bacterial Genomic DNA Extraction Kit (Sangon Biotech). The concentration of the DNA samples was detected using a Qubit fluorescence quantifier; the integrity of the DNA samples was detected by 1% agarose gel electrophoresis. The DNA samples were ultrasonically interrupted using a Covaris instrument, and fragment selection was carried out with an Agencourt AMPure XP-Medium kit to make the sample bands concentrated around 300–400 bp; the amount of purified DNA samples was detected using the kit Qubit dsDNA HS Assay Kit 500 assays. Double-stranded DNA end repair was then performed and an A base was added to the 3’ end. The DNA was ligated to the junction and the ligation product was amplified. The amplification products were screened for fragments with the reagent Agencourt AMPure XP-Medium. PCR products were detected with an Agilent 2100 Bioanalyzer. After denaturing the PCR products to single-stranded forms, single stranded cyclization of the single strands was performed and then the uncyclized linear DNA molecules were digested to obtain the final library. The fragment size and concentration of the library were detected using an Agilent 2100 Bioanalyzer (Agilent DNA 1000 Reagents). The single-stranded circular DNA molecules were replicated by ring rolling to form a DNA nanoball (DNB) containing more than 300 copies. The resulting DNBs were spiked into mesh pores on the chip using high-density DNA nano-chip technology. Sequencing was performed by co-probe anchored polymerization (cPAS). The B. fragilis strain S1-BF-GZ0804 and S2-BF-GZ0302 genomes were sequenced using a Pacbio sequelII and DNBSEQ platform at the Beijing Genomics Institute (BGI, Shenzhen, China). Four SMRT cells Zero-Mode Waveguide arrays of sequencing were used by the PacBio platform to generate the subreads set. PacBio subreads of length < 1 kb were removed. The program Canu was then used for self-correction (Parameters: estn = 24, npruseGrid = 0, corOvlMemory = 4). Draft genomic unitigs, which are uncontested groups of fragments, were assembled using the Canu a high quality corrected circular consensus sequence subreads set. GATK (https://www.broadinstitute.org/gatk/) was used to make single-base corrections to improve the accuracy of the genome sequences (Parameters: -cluster 2 -window 5 -stand_call_conf 50 -stand_emit_conf 10.0 -dcov 200MQ0 > = 4).

Genome component prediction

Gene prediction was performed on the S1-BF-GZ0804 and S2-BF-GZ0302 genome assembly by utilizing glimmer3 (http://www.cbcb.umd.edu/software/glimmer/) with Hidden Markov models. tRNA, rRNA and sRNAs recognition made use of tRNAscan-SE (Parameters:–Spec_ tag(BAOG)–o *. tRNA–f *. tRNA. structure) [25], RNAmmer (Parameters:–s Species–m Type–gff *. rRNA. gff–f *.rRNA. fq), and the Rfam (Parameters:–p blastn–W 7–e 1–v 10000–b 10000–m 8–i subfile–o *. blast. m8) database. The tandem repeats annotation was obtained using the Tandem Repeat Finder (Parameters: 2 7 7 80 10 50 2000 -d -h) (http://tandem.bu.edu/trf/trf.html) and the minisatellite DNA and microsatellite DNA were selected based on the number and length of repeat units. Then the Genomic Island Suite of Tools (GIST) was used for genomicis lands analysis(http://www5.esu.edu/cpsc/bioinfo/software/GIST/) employing the IslandPath-DIOMB, SIGI-HMM, and IslandPicker methods. Insertion Sequence (IS) components were analyzed using ISfinder (database update: 2023-11-10) [26]. Finally, prophage regions were predicted using the PHAge Search Tool (PHAST) web server (http://phast.wishartlab.com/) and CRISPR identification using CRISPRFinder.

Gene annotation and protein classification

The best Blast hit was abstracted using the Blast alignment tool for function annotation. Seven databases including KEGG (Kyoto Encyclopedia of Genes and Genomes), COG (Clusters of Orthologous Groups), NR (Non-Redundant Protein Database databases), Swiss-Prot [27], GO (Gene Ontology), TrEMBL, and EggNOG were used for general function annotation. Four databases for the detection of pathogenicity and drug resistance analysis were also utilized. Virulence factors and resistance genes were identified based on the core dataset in VFDB (Virulence Factors of Pathogenic Bacteria), CARD (The Comprehensive Antibiotic Resistance Database) database, as well as CAZy (Carbohydrate-Active enZYmes Database). Type III secretion system effector proteins were detected using EffectiveT3.

Comparative genomics and phylogenetic analysis

The analysis of nucleic acid synteny was conducted utilizing MUMmer and BLAST to examine the core and pan genes of S1-BF-GZ0804, S2-BF-GZ0302, and a reference strain. The genome sequences of the target bacteria were sequenced in accordance with the genome sequences of the reference bacteria, based on the comparative results obtained from MUMmer (parameter: -b 200 -c 65 --extend -l 20). Subsequently, the protein sets of the target and reference bacteria were compared using BLASTp, with the most favorable comparison result selected for each protein in the library to establish protein pairs (besthit). Ultimately, only those protein pairs that were consistent across both comparisons were retained, with the consistency value for each pair calculated as the average of the consistency values from the two comparisons. Protein gene sets of target and reference strains were clustered by the CD-HIT (parameter: -c 0.5 -n 3 -p 1 -g 1 -d 0) rapid clustering of similar proteins software. The threshold of pairwise identity was set at 50% with a 0.7 length difference cutoff in amino acid. The final gene set of the clustering was used as the Pan gene set. The sequences identified in each sample within the cluster, as per the extracted clustering results, constituted the Core gene set. The gene set unique to each sample was referred to as the Specific gene set. Furthermore, the gene set obtained by excluding both the Core and Specific gene sets from the Pan gene set was termed the Dispensable gene set. The gene family was established using target genes and reference strains, employing multiple software tools in the following manner: protein sequence alignment was performed using BLAST, redundancy was removed through the Solar tool, and gene family clustering was conducted on the alignment results utilizing Hcluster_sg software. The phylogenetic tree was constructed by the TreeBeSTutilizing the NJ method (parameter: treebest nj -b 1000).

Results

Morphological, physiological and biochemical identification

Through isolation and purification, we obtained two strains of B. fragilis from tumor and paracancerous tissues, namely ZY0302 and ZY0804. Both strains were considered Gram-negative as they stained red after Gram staining. When viewed under an optical microscope (100× oil microscope), individual bacteria were visualized as short rods with blunt round ends that did not form spores (Fig. 1A and D). After being cultured on the CDC anaerobic blood plate for 48 h, the two strains formed milky white colonies with similar shapes and sizes (Fig. 1B and E). On the BBE blood plates, the colonies were slightly convex, smooth, and neat. The ZY0302 colonies were slightly larger and the color was brown, while the ZY0804 colonies were grayish white. They both could decompose bile aesculin, turning the culture medium black, a black halo around the colonies was also present (Fig. 1C and F).

Fig. 1

Morphological characteristics of strains ZY0302 and ZY0804. (A and D) After Gram staining, observe the morphology of strains ZY0302 and ZY0804 under an optical microscope (100x) The bar represents a unit length of 10 μm. (B and E) The growth of this strain after 48 h at 37 °C on CDC anaerobic blood agar medium. The bar represents a unit length of 1 cm. (C and F) The growth of this strain after 48 h at 37 °C on BBE blood agar medium. The bar represents a unit length of 1 cm

Physiological and biochemical tests indicated that the two strains could use sucrose, glucose, fructose, maltose, lactose, galactose, mannose, D-ribose, honey disaccharide, cellobiose and aesculin as carbon sources. The strains tested positive in the catalase test, and negative in the nitrate (reduction), gelatin test and nitrate (gas production) test (Table 2). In accordance with Berger’s Manual of Bacterial Identification (eighth edition) [28], the physiological and biochemical characteristics of strains ZY0302 and ZY0804 were consistent with B. fragilis characteristics.

Table 2 Physiological and biochemical properties of B. fragilis ZY0302 and ZY0804

Table 3 Repeat sequence and prophage statistics

Phylogenetic analysis

The 16Sr RNAs of the isolates were sequenced and their lengths were 1423 and 1433 bases, respectively. Seventeen 16Sr RNA sequences that were closely related to the sequences of the isolates were selected as internal taxa from the results of BLAST homology matching analysis (The Per.Ident of the 16Sr RNA sequences of strains ZY0302 and ZY0804 with these sequences were 93.98–99.51% and 92.97–99.51%, respectively, and both had the highest Per.Ident with the sequence of B. fragilis strain W12014D), Porphyromonas gingivalis (ATCC 33277) as an outgroup. Maximum likelihood method concordance analysis based on 16Sr RNA gene sequence alignment showed that strains ZY0302 and ZY0804 had the highest homology with B. fragilis (Fig. 2A). In summary, strains ZY0302 and ZY0804 were identified to belong to B. fragilis. The PCR identification analysis of BFT and its specific typing primers identified strain ZY0302 as bft-1-producing ETBF and strain ZY0804 as NTBF (Fig. 2B).

Fig. 2

Phylogenetic tree and typing identification of strains ZY0302 and ZY0804. (A) The evolutionary history was inferred by using the Maximum Likelihood method and Kimura 2-parameter model. Bootstrap values were based on 1,000 resamplings. Strain Porphyromonas gingivalis ATCC 33,277 was used as an outgrouping. Evolutionary analyses were conducted in MEGA11. (B) Gel electrophoresis images of PCR products amplified based on BFT and its typing primers. M: 2000 DL Marker; 1–3: BFT primer. 4–6: bft-1 primer. 7–9: bft-2 primer. 10–12: bft-3 primer. The products of 1,4,7,10 lanes were produced using DNA from strain ZY0302 as a template, 2,5,8,11 using DNA from strain ZY0804 as a template, and 3,6,9,12 using DEPC water as a template. Imaging was performed using a BIO-RAD gel imaging system

Basic genomic features

ZY0302 (S2-BF-GZ-0302) and ZY0804 (S2-BF-GZ-0804) strains belonging to B. fragilis were isolated from tumor and paracancerous tissues of CRC patients, respectively. The chromosome of each strain was a single circular DNA with sizes of 5,141,346 bp and 5,174,916 bp, with an average G + C content of 43.43% and 43.19%, respectively (Fig. 3), these features were consistent with those of B. fragilis (41–44%) [29]. Both strains contained 18 ribosomal operons (5–23 S-tRNA^Ala-tRNA^Ile-16 S), and all amino acid specificity was provided by 73 tRNA genes. A total of 4,370 genes were predicted in the genome of strain ZY0302 and had an average gene length of 1,041.40 bp, which covered 88.52% of the whole genome sequence. Moreover, 4,369 genes were predicted in the genome of ZY0804, average gene length of 1,048.02 bp, and covered 88.48% of the whole genome sequence. Assessment of the GC skew results revealed that genomic rearrangements and horizontal gene acquisition may have occurred in the ZY0302 and ZY0804 genomes.

Fig. 3

Circular maps of the chromosome of B. fragilis strain ZY0302 and ZY0804. Each circle represents (From outer to inner): 1. Genome Size; 2–3. Forward strand gene and reverse strand gene, colored according to cluster of orthologous groups (COG) classification; 4–5. Forward and Reverse Strand ncRNA, including rRNA (Red), sRNA (Green), tRNA (Yellow); 6. Repeat; 7. G + C content; 8. GC skew (G-C/G + C)

The two types of repetitive elements, DNA repeats and prophages, were present in the genomes of ZY0302 and ZY0804.The DNA repeats included tandem repeats (TR) (Table 3A) and regularly interspaced short palindromic repeats (CRISPR) (Table 3B). 9 and 12 prophages were predicted in both ZY0302 and ZY0804 strains, which accounted for about 4.87% and 6.14% of the total genome, respectively (Table 3C). Insertion sequences (Supplementary Table S1) and transposable elements (8 in ZY0302 and 1 in ZY0804) in the genomes of the two strains were found. However, the presence of plasmids were not detected in the complete genomic sequences of either strain.

DNA restriction and modification (R-M) system

The R-M system is regarded as ubiquitous type of defense mechanism that exists to protect the bacterial cell from invading foreign DNA, especially phage infection [30]. The analysis showed that the ZY0302 genome specified one type I and two type III R-M systems, while the ZY0804 genome specified three type I and one type III R-M systems (Supplementary Table S2). The type I R-M system consisted of three subunits including, HsdR, HsdM, and HsdS, where HsdS determined the DNA-binding specificity of methyltransferases and nucleic acid endonucleases [31]. The ZY0804 genome contained two loci (ZY0804 1855–1858; ZY0804 4192–4195) that were similar to a locus in the NCTC 9343 (BF1839-1842) genome. In the NCTC 9343 genome, this locus showed temporally variable expression of HsdS subunits with eight different possible polypeptide combinations. This in turn provided eight different DNA recognition specificities (Shufflon BB) [31]. In addition, the ZY0302 genome encodes the 5-methylcytosine-specific restriction endonuclease McrBC regulatory subunit McrC. McrBC consists of two polypeptides McrB and McrC. McrB is responsible for GTP binding and hydrolysis as well as DNA binding, while the other McrC is responsible for DNA cleavage. McrBC cleaves specific DNA sequences that have been modified with 5-methylcytosine, hence affording protection to the bacterial cell from foreign DNA [32].

Despite the presence of the R-M system in bacteria, phages have evolved a number of mechanisms to circumvent the limitations of the R-M system, including encoding restriction-resistant proteins and limiting the number of recognition sites in their genomes [33]. From our findings the presence of the gene encoding the restriction endonuclease-resistant protein ArdA in the phage gene of ZY0302 (ZY0302 1853) and the presence of the gene encoding the restriction endonuclease-resistant protein ArdC in the phage gene of ZY0804 (ZY0804 2152) may help elucidate the reason why these prophages can successfully be transferred into the bacterial genome.

Comparative genomics

The predicted proteome levels of ZY0302 and ZY0804 were compared with those of B. fragilis obtained through the Genebank and a phylogenetic evolutionary tree was then constructed. ZY0302 and ZY0804 shared 2,775 core genes with the reference strain and possessed 320 and 285 unique genes, respectively (Fig. 4A). Phylogenetic analysis of the core genes revealed that ZY0302 and ZY0804 were in the same evolutionary branch (Fig. 4B). However, the dispensable genes showed similar gene heat in different branches, which was unlike that of the reference strain (Fig. 4C).

Fig. 4

Common and unique genes of ZY0302 and ZY0804 with reference strains. (A) Venn diagram of the Pan gene set. (B) Phylogenetic trees were constructed based on core genes of the sample and reference strains. The NJ algorithm of TreeBeST was used to construct the phylogenetic tree. (C) Dispensable gene heatmap

When compared to the reference strain, gene families showed that in the ZY0302 predicted genome, 4,122 (94.32%) genes were clustered into 2,623 gene families, of which 17 unique paralogous genes belonged to 8 strain unique gene families. The 2,615 immediate homologous gene families included 1,962 multicopy immediate homologous genes. In the ZY0804 genome, 4,152 (95.03%) genes were clustered into 2,694 gene families, with 18 unique paralogous homologs in 9 strain unique gene families and 1,948 multicopy orthologs in 2,685 orthologous gene families. The largest lineal homologous family was composed of them while the reference strain consisted of 471 homologous genes, including genes such as BF9343_RS01055 and BF_RS10050, which are RagB/SusD family vegetative uptake outer membrane protein genes. Nucleic acid collinearity results showed that the genomes between strains ZY0302 and ZY0804 had similar structural variability during the evolutionary process, with different degrees of genomic rearrangements, inversions, insertions, and deletions with different reference strains (Supplementary Fig. S1).

Carbohydrate metabolism

B. fragilis can utilize host-derived cell surface glycoproteins and glycolipids as an energy source, and these may include monosaccharides (e.g., glucose, galactose, mannose, etc.) and more complex compounds (for example, N-acetyl-d-glucosamine and N-acetylneuraminic acid.) [6]. Annotation using the CAZy database revealed the presence of various carbohydrate-associated enzymes in the genomes of ZY0302 and ZY0804, including Glycoside Hydrolases(GHs), GlycosylTransferases (GTs), Polysaccharide Lyases (PLs), Carbohydrate Esterases(CEs), Auxiliary Activities(AAs), and Carbohydrate-Binding Modules(CBMs) (Table 4).

Table 4 Taxonomic annotation of carbohydrate enzymes of ZY0302 and ZY0804

Virulence factors

By comparing the ZY0302 and ZY0804 genomes with the Virulence Factor Database (VFDB), 123 and 121 virulence genes were predicted, respectively (Supplementary Table S3). Consistent with the PCR results, the 16 virulence factors specific to ZY0302 included zinc-dependent metalloproteinase toxin (bft1), the major virulence factor of ETBF. Among the 11 specific virulence factors of ZY0804, α-1, 3-fucosyltransferase was identified, this enzyme is a key to the last step of the formation of LPS Lewis antigen. Fucose-containing LPS are commonly expressed by pathogenic bacteria such as Helicobacter pylori and E. coli, and they are thought to be involved in molecular mimicry, adhesion, colonization, and modulation of the host immune response [34].

Although B. fragilis is exclusively anaerobic, it is also highly tolerant to oxygen and can survive exposure to an aerobic environment for 48–72 h [35]. The high tolerance of B. fragilis to oxidative stress is thought to be an important factor in its escape from the anaerobic environment of the gastrointestinal tract to aerobic tissues such as the peritoneal cavity as well as its tolerance to the respiratory burst induced by immune cells [35]. We found the presence of catalase, superoxide dismutase, aspartate 1-decarboxylase and peptide methionine sulfoxide reductase msrA/ msb in the genomes of ZY0302 and ZY0804, as well as alkyl hydroperoxide reductase AHPC and the heat shock protein CLPB. Among them, CLPB is a member of the AAA family (ATPases associated with a wide range of cellular activities), which has been shown to be important in a variety of pathogenic bacteria and is an important virulence factor, playing a key role in the cellular response to stressful conditions such as heat stress, antibiotic resistance, and oxidative stress [36]. In addition, CLPB proteins regulate the secretion of bacterial effector molecules associated with the type VI secretion system (T6SS), which is critical for the survival and infectivity of clinically relevant microorganisms [36]. However, there are no studies that directly refer to the specific functions or roles of CLPB proteins in B. fragilis.

Among the ZY0302 and ZY0804 virulence genes, we identified genes related to type IV bacteriophage biosynthesis, the hmu direct heme uptake system (two in ZY0302 and one in ZY0804) and PDTC (pyridine-2,6-dithiocarboxylic acid). PDTC is thought to act as an iron carrier in Pseudomonas aeruginosa and is associated with the transport of iron and other metals [37]. Iron is removed from heme in the periplasmic space and the released iron needs to be transported to the cytoplasm of B. fragilis via the FeoAB system and the CobN-like proteins BtuS1 and BtuS2 [38]. It is therefore believed that B. fragilis employs similar systems in its virulence mechanism. Through COG database gene function annotation, we found the presence of FEOB (ZY0302 2050; ZY0804 1995) and COBN proteins (ZY0302 3039, 3304; ZY0804 3338) in the genomes of ZY0302 and ZY0804. In addition, five VirB4 subunits of the type IV secretion system (1867, 1969, 3147, 4060, 4109), six VirD4 subunits (1054, 1876, 1979, 3134, 4023, 4071), and three TrbF fractions were identified in the genome of ZY0302 by functional annotation of genes in the COG database (1864, 1966, 4056), and the genome of ZY0804 identified one VirB4 subunit (2067) and one VirD4 subunit (2085) of the type IV secretion system. Further analysis showed that most of these subunits were not present on the bacterial genome chromosome but originated from prophages, except for one VirB4 subunit (3147) and three VirD4 subunits (1054, 3134, 4023) which were identified on the bacterial chromosome of the ZY0302 genome.

B. fragilis has been shown to induce abscess formation and its potential to induce abscess formation is significantly related to its capsule [39]. Virulence gene analysis showed that there were 74 (74/121) virulence genes were involved in capsule formation in ZY0804 and 76 (76/123) in ZY0302. At least eight and nine capsular biosynthesis gene clusters in their genomes were found in the ZY0804 (107 to 112meme 1442 to 1447 meme 1456 to 1465, 2038 to 2046, 2253 to 22268, 3228 to 3232, 3516 to 3518, 4170 to 4172) and ZY0302 (95 to 99991436 to 1447, 1701 to 1701 to 1776 to 2087, 2183 to 2185, 2499 to 2502, 3192 to 3194, 3488 to 3488 to 4215), respectively. ZY0804 had two capsular biosynthetic gene clusters (2038 to 2046 and 3228 to 3232) derived from prophages and ZY0302 had one (2076 to 2087).

B. fragilis can produce a variety of extracellular and physiological enzymes to attack and penetrate the host’s extracellular matrix. Extracellular enzyme-related genes such as neuraminidase, lipase, hemolysin (Table 5), protease, and phosphatase (not listed) were predicted in the genomes of ZY0302 and ZY0804 using COG database gene function annotation. The products synthesized by these extracellular enzyme genes play important roles in the adhesion, invasion and growth of B. fragilis.

Table 5 Extracellular enzymes of ZY0302 and ZY0804

Discussion

In this study, we isolated ZY0302 and ZY0804 from tumor tissues and paracancerous tissues of CRC patients. They were confirmed to be B. fragilis by morphological analysis, physiological and biochemical characteristics, and 16 S rRNA sequencing according to Berger’s Manual of Bacterial Identification. Furthermore, we detected the expression of BFT gene in the two bacterial strains by applying PCR. The final strains, ZY0302, were identified as a bft-1-producing ETBF, and ZY0804 as an NTBF. Whole genome sequencing (WGS) is able to identify taxa at a higher taxonomic resolution relative to 16 S rRNA and provides more information about gene function. We characterized the genomes of the isolates and identified and analyzed their putative virulence factors.

Pan-genomic analysis revealed that the isolates shared multiple immediate homologous gene families with the B. fragilis reference strain. RagB/SusD family nutrition uptake outer membrane protein gene was observed to be their largest gene family. SusD-like proteins have α-helix folding, bind polysaccharides on the cell surface, and may promote polysaccharides to shuttle to related TonB-dependent porins for transport to periplasm for degradation [40]. Also, the SusC family of outer membrane proteins is the largest paralogous family of B. fragilis and Bacteroides pleomorphic [19]. The multiple copies of these polysaccharide utilization genes, as well as the presence of large numbers of GHs, GTs, and CBMs in the genome, may have been designed to assist in the effective utilization of nutrients in the distal gut. The distal gut is deficient in monosaccharides and disaccharides, and these amenities could help in the efficient utilization of polysaccharides and a variety of dietary carbohydrates. Phylogenetic analysis of core and dispensable genome sequences showed that B. fragilis isolated from cancerous and paracancerous tissues belonged to the same evolutionary branch but was in a different evolutionary branch from B. fragilis from other tissue sources (BF-BE1 was isolated from a patient with a wound infection, BF-BOB25 was isolated from a fecal sample from a patient with intestinal ecological dysbiosis, BF-NCTC − 9343 was isolated from patients with abdominal infections, BF-Q1F2 was isolated from patients with skin and soft tissue infections, and BF-YCH46 was isolated from patients with sepsis). From these results, we speculated that ZY0302 and ZY0804 co-evolved the core and non-essential genomes as a means to adapt to the changing microenvironment of the tumor and paracancerous tissues, thus enabling them to colonize the ever-changing microenvironment long-term. In addition, ZY0302 and ZY0804 have undergone structural variation in gene rearrangements (inversions, insertions, and deletions) during evolution. This crucial evolutionary driver regulates gene expression and advances ecological niche adaptation [41]. Through a process known as horizontal gene transfer (HGT), some virulence factors and antibiotic resistance genes that were not originally present in the bacterial genome may be acquired. Consequently, horizontally transferred genes in a recipient bacterium may incite a new adaptation thereby affording the bacterium new ecological niches that were inaccessible through mutation alone [42]. The presence of phages, a class of mobile genetic elements that can be transmitted from cell to cell, may also allow some bacteria to acquire antibiotic resistance, enhance adaptation to the environment, improve adhesion and enable bacteria to become pathogenic [43]. Insertion Sequences (ISs), transposases, and mobile elements from phages were identified within the genomes of ZY0302 and ZY0804. The phages contained several genes associated with capsular biosynthesis and the type IV secretion system, while they lacked other virulence genes and antibiotic resistance genes. The presence of these phages enables better adaptation of the strain to new ecological niches and colonization.

B. fragilis makes up the gut microbiota and remains commensal with the host gut, however under unfavorable conditions, this bacteria can turn into an opportunistic pathogen leading to abscess formation and bacteremia in several parts of the body (abdomen, brain, liver, pelvis, and lungs) [5]. Several virulence factors contribute to its conversion from a symbiont into a pathogenic bacterium, including complex capsular polysaccharides, BFT, LPS, and virulence-related enzymes [6]. Furthermore, recent studies have shown that the B. fragilis correlation with the development of CRC [18, 44]. The development of CRC is accompanied by the generation of an inflammatory microenvironment and excessive oxidative stress [45]. ZY0302 and ZY0804 encode multiple genes involved in the regulation of oxidative stress, which is favorable for their survival in the excessively oxidatively stressed tumor microenvironment.

B. fragilis effectively evades the host’s immune response due to its unique cellular structure. It possesses the ability to modulate the surface composition of its capsular polysaccharides, which contributes to its capsule’s role in increasing the bacteria’s resistance to complement-mediated lysis as well as to phagocytosis and subsequent destruction [19, 20, 46, 47]. In addition, its capsule and LPS also act as adhesins, allowing the bacteria to attach and colonize the infection site [6]. ZY0302 and ZY0804 from CRC and paracancerous tissues are rich in capsular biosynthesis genes which is consistent with the characteristics of B. fragilis. Among the predicted virulence factors, 61.79% of them for ZY0302 and 61.16% for ZY0804 are related to capsular biosynthesis, and have multiple capsular biosynthesis gene clusters. Certain gene clusters originate from phages via horizontal gene transfer, encoding proteins associated with cell structure derived from phage sequences. This process can further augment the adaptability of bacteria to their environment and enhance their adhesion capabilities. We also identified some genes in the genomes of ZY0302 and ZY0804 that are related to the synthesis of type IV bacterial hairs. B. fragilis can use bacterial hairs to adhere to other bacterial species and host tissues, which is important for bacterial mucosal colonization [48]. In addition, the genomes of ZY0302 and ZY0804 contain the kdsA gene encoding 2-dehydro3-deoxyphosphate octanoate aldolase. KdsA genes have been implicated in LPS synthesis with adhesion and endotoxin properties [49]. B. fragilis strains isolated from patients with polyps are enriched in the LPS gene and have a greater ability to activate TLR4 and induce a proinflammatory cytokine response [18]. ZY0302 and ZY0804 encode genes for glycosyltransferases involved in LPS biosynthesis, and 18 and 16 of their virulence genes, respectively, are associated with LPS synthesis (Supplementary Table S3).

Iron and heme are essential nutrients for the growth and survival of Bacteroides, and key determinants of bacterial virulence as well. The hmu loci identified in the ZY0302 and ZY0804 genomes play a major role in heme acquisition in Porphyromonas gingivalis, and a similar locus was found in B. fragilis NCTC9343 [50]. Alluding to its vital potential role in heme acquisition of B. fragilis. Bacterial type IV secretion systems (T4SSs) are a functionally diverse translocation superfamily that not only mediate the horizontal transfer of DNA, but also transport macromolecules such as proteins and ribonucleoprotein complexes, which allows the bacteria to encode resistance to heavy metals and antibiotics [51]. However, previous studies have shown that there are no traces of III type, IV type, autotransporter or two-partner secretion systems in the genome of B. fragilis [52]. In this study, we discovered that the loci linked to the type IV secretion system were identified within the genomes of ZY0302 and ZY0804; however, they were predominantly located within the phage genes, suggesting that these genes were likely obtained through horizontal gene transfer.

B. fragilis can produce a variety of extracellular and physiological enzymes, with neuraminidase being the most prevalent. This enzyme cleaves mucin polysaccharides, thereby promoting microbial growth, and serves as a significant virulence factor that impairs the immune function of the host [6]. Lipase is a secreted virulence factor, and the putative roles of microbial secretion of lipase during infection include host cell adhesion, lipid digestion for nutrients and triggering of the inflammatory cascade [53]. The production of hemolysins may facilitate the in vivo access of B. fragilis to iron and heme through the destruction of host cells and erythrocytes. Two hemolysins, namely HlyA and HlyB, have been identified in B. fragilis, exhibiting synergistic functions in the process of erythrocyte hemolysis [54]. The proteases attack a range of host proteins, including extracellular matrix proteins and cell adhesion molecules, and their disruption leads to the loss of cell surface receptors and tissue integrity [55, 56]. It is evident that these enzymes not only play a key role in their pathogenicity, but are also involved in bacterial nutrient metabolism and immune escape to enhance their own survival.

The genomic diversity and genetic differences between strains of B. fragilis may be related to its adaptations in different hosts and environments. Considering the dual functions of B. fragilis in both intestinal health and disease, it is imperative that future research delves deeper into the ecological implications of these genomic variations and their impact on host health and disease. Utilizing high-throughput sequencing technology, researchers should perform extensive genomic comparisons of B. fragilis strains derived from various sources to uncover their evolutionary connections and functional distinctions. This approach will lay the groundwork for the formulation of targeted intervention strategies. However, the limited sample size of B. fragilis examined in this study may have led to findings that do not provide robust evidence regarding their epidemiological significance. The primary factor contributing to this limitation is that, in the present study, we isolated and identified B. fragilis from a single tissue section. This approach effectively addresses issues of contamination and pollution, thereby ensuring the specificity of the strain source. However, it also heightens the risk of leakage, which may be a significant reason for our limited positive sample size. Future studies could enhance this sample size by either homogenizing the tissue or developing a more systematic protocol to improve the isolation efficiency of B. fragilis. In addition, although we observed the correlation between genomic alterations in B. fragilis and adaptation to the tumor microenvironment, we did not directly confirm that these alterations promote B. fragilis adaptation to the tumor microenvironment through corresponding experiments. In future studies, we will also further validate the role and significance of these genomic alterations in the adaptation of B. fragilis to the tumor microenvironment in combination with modern genetic engineering techniques such as CRISPR-Cas9.

Conclusions

The genomic genetic characteristics of B. fragilis strains ZY0302 and ZY0804 isolated from tumor and paracancerous tissues of CRC patients are different from those isolated from non-tumor tissues, the selection pressure from the tumor microenvironment may be the important driver of pan-genomic variability in B. fragilis. Changes in the genomic characteristics of ZY0302 and ZY0804 may facilitate bacterial adaptation to the changing tumor microenvironment, where gene rearrangement and horizontal gene transfer offer the possibility of better adaptation to changes in the microenvironment. Of course, These alterations may also elevate the infectious potential of the bacteria; however, further research is necessary to elucidate their roles in the progression of CRC.