Bacteroides uniformis CECT 7771 was isolated from stools of a healthy infant included in a prospective observational study carried out in a sub-group of 75 full-term newborns with at least one first-degree relative with celiac disease belonging to the Proficel study cohort (Sanchez et al., 2011). Infants were enrolled in the study after written informed consent was obtained from their parents as indicated elsewhere (Sanchez et al., 2011). B. uniformis CECT 7771 was grown for 48 h in Schaedler anaerobe agar plates (Oxoid) supplemented with 100 mg/L kanamycin, 7.5 mg/L vancomycin, and 0.5 mg/L vitamin K₃. Bacteria were incubated at 37°C under anaerobic conditions into a Whitley DG250 Anaerobic Workstation (don Whitley Scientific, Inc., Shipley, UK). Single colonies were used to inoculate 10 mL of modified Schaedler broth (10 g/L tryptone soy broth, 2.43 g/L casein pancreatic peptone, 0.43 g/L soy peptone, 2.15 g/L meat extract, 5 g/L yeast extract, 5 g/L glucose, 0.75 g/L Tris-HCl, 0.4 g/L L-cysteine, 0.01 g/L hemin and 0.5 mg/L vitamin k₃). The L-cysteine, hemin and vitamin K were filtered individually by using 0.22 μm disposable filters (Millipore) and then were added to the autoclaved media. Ten mL of an overnight culture of the strain studied was used for genomic DNA isolation. To evaluate the effects of different carbon sources on B. uniformis CECT 7771 growth, overnight cultures (0.7–0.8 OD₆₀₀) of this bacterial strain were diluted 1/20 in fresh pre-warmed and oxygen-depleted modified Schadler media containing 0.5% w/v of the different carbon sources including glucose (Scharlau, Cat#GL01271000), inulin (Sigma, Cat#I2255), wheat bran extract (WBE) (Cargill, Antwerp, Belgium), gum arabic (Sigma, Cat#30888), pectin (Sigma, Cat#P9135), or type II mucin from porcine stomach (Sigma, Cat#M2378) and growth kinetics were monitored for 10 h. Each carbon source stock was sterilized by filtration using 0.22 μm disposable filters, thus avoiding the alteration of their original structure and chemical composition. The OD₆₀₀ was measured at 1-h intervals for each sample until stationary phase. An aliquot of respective cultures was collected at exponential growth phase (OD₆₀₀ ~0.4) for RNAseq analysis. The cells were then pelleted by centrifugation at 3,000 × g for 15 min at 4°C. The supernatant was aspired off, filtered using 0.22 μ filters and immediately stored at −80°C for further metabolite analyses and the cell pellet was stored at −80°C for RNA extraction. Three independent experimental trials were included in the analyses.

DNA and RNA from cell cultures were isolated using the MasterPure^TMGram Positive DNA Purification Kit (Epicenter) with slight variations over manufacturer's instructions. Briefly, a cell lysis step was improved by incubating each cell suspension with 500 mg Lysozyme (Sigma, Cat #62970) and 20 U Mutanolysin (Sigma, Cat #9901) for 60 min at 37°C. For RNA isolation samples were incubated with 2 U DNase I (Epicenter) at 37°C for 60 min instead of the RNase A treatment.

B. uniformis CECT 7771 DNA was sequenced using MinION portable DNA sequencer and flowcells based on the R7.3 and R9.4 pore chemistries. The R7.3 run consisted of a DNA library prepared with the Genomic DNA Sequencing Kit SQK-MAP006 (Oxford Nanopore Technologies) and 6 μg genomic DNA according to manufacturer's instructions. The genomic DNA was sheared with gTUBE^TM (Covaris) to get fragments of about 10 Kbp, the sheared DNA was then repaired with the PreCR® Repair Mix (New England Biolabs) followed by purification using Agencourt AMPure XP beads (Beckman Coulter). DNA attached to AMPure beads was washed twice with freshly prepared 70% ethanol and magnetic rack. The DNA was eluted in 46 μL nuclease-free water and quantified by using Qubit 3.0 fluorometer and the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). One microgram of repaired DNA was further processed using the NEBNext Ultra II End Repair/dA-Tailing Module (New England Biolabs) followed by a new round of washing in the magnetic rack as previously done. The DNA was eluted from magnetic beads using 38 μL nuclease-free water and adapter ligation step was conducted by adding 10 uL Adapter Mix, 2 uL HP adapter, and 50 uL Blunt/TA Master Mix (New England Biolabs) mixing by inversion between each sequential addition. After 15 min incubation at room temperature 1 μL HP Tether was added to ligation reaction and incubation was extended for 10 min. The adapter-ligated genomic DNA was recovered with Dynabeads® MyOne Streptavidin C1 beads (Thermo Fisher Scientific) a magnetic rack, and Elution Buffer (SQK-MAP006 kit). Approximately, 205 ng DNA library were recovered and they were loaded into a brand new, sealed R7.3 flowcell previously fitted to the MinION™and primed twice with 71 μL premixed nuclease free water, 75 μL 2x running buffer, and 4 μL fuel mix. The initial sequencing mix was prepared with 56 μL nuclease-free water, 75 μL 2X running buffer, 4 μL fuel mix, and 15 μL DNA library (~120 ng). A standard 48-h sequencing protocol was initiated using MinKNOW^TM v0.50.2.15. The base-calling was performed through data transference using the Metrichor™ agent v2.36.2 and 2D Basecalling workflow v1.62. During the sequencing run, the remained aliquot of DNA library was loaded in equally conditions after 24 h of initial input. The R9.4 run consisted of a DNA library prepared with the 1D Genomic DNA Sequencing Kit SQK-LSK108 (Oxford Nanopore Technologies) and 3 μg genomic DNA according to manufacturer's instructions. Briefly, the DNA repair kit used was NEBNext FFPE RepairMix (New England Biolabs) and the final DNA library purification was performed with AMPure beads instead Dynabeads®. The sequencing mix was obtained by combining 37.5 μL running buffer, 25.5 μL LLB beads, and 12 μL DNA library (~240 ng) and loaded into the spot-on port of R9.4 flowcell previously mounted into a MinION^TM MkIb sequencer. A standard 48-h sequencing protocol was initiated using MinKNOW^TM v1.1.17. The base-calling was performed through data transference using the Metrichor™agent v2.43.1 and 1D BasecallingFLO-MIN106 450bs workflow v1.121. During the sequencing run, the remained aliquot of DNA library was loaded in equally conditions after 12 h of initial input. This sequencing run was extended for 24 h in sum.

The total RNA obtained was quantified by using the Qubit 3.0 fluorometer and Qubit® dsDNA HS Assay Kit (Thermo Fisher Scientific). Thus, 30 μg total RNA were obtained by pooling respective replicates of different culture conditions in equimolar quantities (10 μg per replicate). RNA pools were sent to Eurofins Genomics GmbH (Ebersberg, Germany) to produce cDNA libraries with an insert size of ~400 bp with prior rRNA depletion using RiboZeroTM Magnetic Kit Gram-Negative Bacteria (Epicenter). The six cDNA libraries (glucose, WBE, pectin, inulin, gum arabic, and mucin) were pooled and sequenced in one Illumina HiSeq2500 channel with chemistry v4 and configuration 2 × 125 paired-end reads.

Quality assessment of fast5 files and conversion to fastq and fasta formats was performed using the poRe (Watson et al., 2014) package v0.17. The B. uniformis CECT 7771 genome assembly was performed with 2D reads from R7.3 chemistry (32,764 reads), the high-quality set of 1D reads from R9.4 chemistry (10,598 reads), and the Canu assembler v1.3 (Koren et al., 2017) with options genomeSize = 4.5 m, minReadLength = 500, corMinCoverage = 5. A genome sequence refinement was performed with the RNAseq data derived from Illumina HiSeq2500 consisting of ~275 million paired-end reads, which were mapped against the draft genome of B. uniformis CECT 7771 using bowtie2 (Langmead and Salzberg, 2012) for aligning and samtools v1.3.1 (Li et al., 2009) set of algorithms for indexing, sorting, and pileup of mapped reads. A consensus sequence was recovered and further processed to allow super-scaffolding. This last step was completed by performing de novo transcriptome assembly using the entire set of RNAseq paired-end reads and the velvet v1.2.10 and oases v0.2.09 assemblers using 63 as k-mer parameter (Zerbino and Birney, 2008; Schulz et al., 2012). Thus, the entire set of contigs assembled (31,576 in total) were mapped against the draft genome of B. uniformis CECT 7771 obtaining a new and more refined consensus sequence. Finally, a blast-based comparison among refined scaffolds of B. uniformis CECT 7771 was carried out to detect strong similarity among edges. This information was used to perform manual assembly of the 11 initial scaffolds generated from Canu into 6 superscaffolds. The genome assembly of the B. uniformis CECT 7771 was submitted to the European Nucleotide Archive (ENA) where it is publicly available under primary accession number PRJEB19372 (ENA WGS accessions: FZQS01000001-FZQS01000006). Comparative genomics against the B. uniformis ATCC 8492 (ENA WGS accession: AAYH02000000), B. uniformis CL03T00C23 (ENA WGS accession: AGXY01000000), B. uniformis dnLKV2 (ENA WGS accession: ASSO01000000), and B. uniformis 3978-T3i (ENA WGS accession: JNHO01000000) was assessed with alignment algorithms implemented in Mauve v2.3.1 (Darling et al., 2010) and BRIG (Alikhan et al., 2011). Furthermore, genome-wide average nucleotide identity approach implemented in JSpecies v1.2.1 (Richter and Rossello-Mora, 2009) was used to reconstruct phylogenetic relationships among B. uniformis species and other Bacteroides such as Bacteroides vulgatus ATCC 8482 (ENA accession: CP000139.1), Bacteroides fragilis NCTC 9343 (ENA accession: CR626927.1), Bacteroides thetaiotaomicron VPI 5482 (ENA accession: AE015928.1), Bacteroides caccae ATCC 43185 (ENA WGS accession: AAVM02000000), and Bacteroides stercoris ATCC 43183 (ENA WGS accession: ABFZ02000000). Values from ANIm and tetranucleotide distribution were used to reconstruct a UPGMA dendrogram (http://genomes.urv.cat/UPGMA) using the RSMD distance coefficient and the iTOL web server (Letunic and Bork, 2016) was employed to draw the respective rooted trees. Gene prediction and functional annotation were achieved by using Prodigal v2.6 (Hyatt et al., 2010), tRNAscan v1.4 (Lowe and Eddy, 1997), RNAmmer v1.2 (Lagesen et al., 2007), CRISPR Recognition Tool (CRT) (Bland et al., 2007), KEGG Automatic Annotation system (Moriya et al., 2007), CAZy database (Lombard et al., 2014), SMART database (Letunic et al., 2014), and CAT server (Park et al., 2010). Venn diagrams were designed in jvenn server (Bardou et al., 2014). Assessment for presence of virulence factors, antibiotic resistance genes, and toxins was done by executing local blast to compare all the ORFs predicted against the annotated proteins present in the VFDB (Chen et al., 2016), MvirDB (Zhou et al., 2007), and ARDB (Liu and Pop, 2009) databases. Processing of RNAseq data was done according to previous studies (Benitez-Paez et al., 2016b). Briefly, the quality filtering and trimming were performed using FASTX-toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). Read mapping was assisted by Blast algorithm (Altschul et al., 1990) and selecting alignments >50% of read length (>70 nt) and 100% identity. Read counts were normalized using RPKM and a genome-wide exploratory differential expression among transcriptomes generated from usage of different carbon sources was measured with GFOLD (Feng et al., 2012) using the glucose-derived expression as baseline feature and followed by further evaluation with qPCR methods. In order to increase the stringency for detecting more probable signals of differential expression, we only selected genes with GFOLD score ≤ −1.5 or ≥1.5 (log₂ fold-change) despite than any gene with GFOLD score different than zero would be indicative of up- or down-regulation. Sequence information supporting the six B. uniformis CECT 7771 transcriptome analyses was submitted to the ENA where it is publicly available under primary accession number PRJEB19372. Hierarchical clustering of genes differentially expressed in different culture conditions was achieved by using Euclidean distance and average linkage methods (Metsalu and Vilo, 2015).

The genes BUNIF7771_0387, BUNIF7771_0544, BUNIF7771_0548, BUNIF7771_1668, BUNIF7771_1883, BUNIF7771_3473, BUNIF7771_3732, and BUNIF7771_4131 were selected from the exploratory RNAseq analysis to assess specific changes in expression by qPCR. The gene-specific oligonucleotides used for this aim are presented in the Table S1. The cDNA was synthesized using 5 μg of total and non-pooled RNA remaining from that used for the RNAseq approach (three replicates per treatment), and the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems) according to the manufacturer's instructions. The qPCR reactions were set in 96-well plates using the SYBR Green I Master Mix (Roche Lifesciences), 0.5 μM of forward oligonucleotide, 0.25 μM of reverse oligonucleotide, and 1 μL of the cDNA reaction. All treatment samples were set in triplicate in the plate and amplified in a LightCycler 480 II with the following cycling profile: initial incubation at 95° for 5 min and 40 cycles of 10 s at 95°, 20 s at 63°, and 15 s at 72°. Finally, the melting curve was set from 65 to 97° with a ramp rate of 0.11°/s. The expression level for each gene was measure according to the ΔΔCt method, using the expression of the 16S rRNA gene as calibrator, and expression of glucose samples as reference. RQ values were finally obtained with calculation of 2^−ΔΔCt for all samples and replicates. Differential expression was assessed by the one-sided t-test with Welch's correction supporting pairwise comparisons between gene expression under glucose and remaining treatments. The copy number of the pBU7771 extrachromosomal element was calculated by absolute quantification using the primers pBU7771-F and pBU7771-R (Table S1). The single-stranded DNA (ssDNA), fully covering the region to be amplified (111 nt) was obtained from Isogen Life Science B.V (Utrecht, The Netherlands) where it was synthesized, PAGE-purified, and quantified, and used in molecule titration during qPCR. The number of plasmid molecules per ng DNA was obtained and divided by the theoretical number of B. uniformis CECT 7771 genomes presented in 1 ng DNA (5.16 Mb) (http://cels.uri.edu/gsc/cndna.html), this ratio was used to infer the number of plasmids per cell.

The gamma-amino butyric acid (GABA) concentration was analyzed in supernatants (150 μL) obtained from cultures supplemented with either glucose, mucin or pectin and collected at OD₆₀₀ ~ 0.4 by LC/MS approach in the Central Service for Experimental Research (SCSIE) at the University of Valencia. The GABA (Sigma, Cat#A2129) standards were prepared in milliQ grade water. Both the standards and samples were diluted 1:4 with acetonitrile in 200 μL reaction volume, then samples were centrifuged at 10,000 × g for 10 min to facilitate protein precipitation and 50 μL of supernatant were loaded in the ACQUITY® TQD (Waters Corporation. Milford, MS, USA) LC instrument coupled to an electrospray ionization (ESI). GABA was separated in a HILIC Kinetex column (Phenomenex®) using 0.1% formic acid: 75% acetonitrile at 23°C. The area under the curve (AUC) of GABA standards at 1, 0.1, 0.01, and 0.001 mM were used to quantify its concentration in the culture supernatants. Statistical analysis was done with the values of three independent assessments using one-way ANOVA with pairwise comparisons and Bonferroni correction.

Following an extensive sequence analysis of the genomic information retrieved from a third generation sequencing platform, the MinION^TM, we could initially assemble the B. uniformis CECT 7771 genome which consisted of ~5.16 Mbp with a N50 of 1.32 Mbp arranged in 11 scaffolds. The MinION^TM data is characterized by producing error-prone DNA reads with a per-base accuracy above 85% in 2D reads (Ip et al., 2015), therefore we assumed further and hybrid approaches to get genome sequence refinement. After refinement steps with RNAseq data, using individual reads and assembled into longer transcript contigs, derived from six different experiments to assess the transcriptional profile under different carbon sources, we obtained a high quality genome assembly consisting of 5.16 Mpb with a N50 2.44 Mpb and arranged in 5 super-scaffolds representing the chromosomal genetic information of B. uniformis CECT 7771. This genome information was submitted to the European Nucleotide Archive (ENA) and can be accessed with identifiers FZQS01000001 to FZQS01000005. Additionally, we found an extra-chromosomal element of 2,746 bp in length that seems to encode neither virulence nor antibiotic resistance genes (ENA accession FZQS01000006). By absolute quantification approaches, we found the pBU7771 plasmid is present in very-low copy number reaching a proportion of 1.86 plasmids per cell. We have predicted that in the current version of the B. uniformis CECT 7771 genome there are 5,226 ORFs, 4 ribosomal RNA operons, and 67 tRNAs encoded. Moreover, we found one CRISPR region characterized to have 11 repeats consisting of the following DNA sequence: GTTGTGATTTGCTTTCATTTTAGTATCTTTGAACCATTGGAAACAGC and 10 spacers with an average length of 30 nt. This last feature was confirmed by the fact that the downstream of this CRISPR region are contiguously located the BUNIF7771_0521 and BUNIF7771_0522 ORFs encoding the Cas2 and Cas1 CRISPR-associated endonucleases, respectively. Ribosomal RNA sequence information of B. uniformis CECT 7771 is publicly available in the ENA under accession numbers LT745888, LT745889, and LT745890 for 5S, 16S, and 23S molecules, respectively. An initial comparison with other genome assemblies from B. uniformis strains is presented in Table 1. The genome assembly presented in this study fully fits the features observed for previously sequenced B. uniformis strains. Notably, the draft genome of B. uniformis CECT 7771 share similar GC content with its counterparts being 46.7% and constitute the largest genome assembly of any B. uniformis published with the lowest number of scaffolds reconstructed. In this regard, it seems that the hybrid assembly strategy based on long read sequencing approach followed by refinement with massive short read data resulted in a valid methodology to sequencing the B. uniformis CECT 7771 genome retrieving more sequence information to distinguish interesting features of B. uniformis strains. Among such features, we highlighted the presence of 4 complete ribosomal RNA operons whose number was not consensual from previous assemblies. Besides, we reported for the first time the presence of a plasmid element in a B. uniformis strain of 2,746 bp in length. Using this element as query in a Blast searching against the NCBI's non-redundant nucleotide database we observed that this element was commonly found in other Bacteroides species consisting of an extrachromosomal DNA element that was also recovered from metamobilome studies in rat cecum (Jorgensen et al., 2014) with a high sequence similarity to a plasmid element borne by B. fragilis IB143 (Smith et al., 1995). Notwithstanding, none of the three ORFs predicted to be encoded by this element match with any antibiotic resistance, virulence, or pathogenic gene annotated in the ARDB, MvirDB, and VFDB databases.

Comparative genomics approaches based on whole-genome alignment were used to identify distinctive features of B. uniformis CECT 7771 compared to other strains belonging to the same and to other species. The average nucleotide identity and tetranucleotide distribution across the genome indicated that B. uniformis CECT appears to be very related to other B. uniformis strains but its genomic structure also differs from that of other Bacteroides species (Figure 1 and Table S2). A more detailed comparison among different B. uniformis strains has led to identify genomic regions exclusively present in the genome of B. uniformis CECT 7771, named specificity regions hereinafter (SR1–SR5, Figure 2). Despite the low level of functional annotation recovered when ORFs encoded into the major SR regions were analyzed by KAAS server (Moriya et al., 2007), we could identify interesting features such as presence of endopeptidases (SR1 and SR3), DNA methylases (SR1, SR2, and SR3), a type-IV restriction endonuclease (SR3), a two-component system (SR1), and polysaccharide metabolism associated enzymes (SR2) for importing, degradation and synthesis.

The members of the phylum Bacteroidetes are well known components of the gut microbiota thriving in environments enriched in glycans. The sus genes, an archetypal locus of membrane proteins for binding and importing polysaccharides in Bacteroidetes species, are a fundamental part of the PULs. The number of these loci specialized in glycan degradation in Bacteroides spp. and Prevotella spp. species outnumbers the number of PULs found in species of other phyla such as Firmicutes (White et al., 2014). Therefore, we wanted to survey the number and classes of enzymes involved in binding, degradation, modification, and synthesis of glycans in the B. uniformis CECT 7771 genome. Using the CAT server (Park et al., 2010), we retrieved annotation of ORFs encoding enzymes homologs to reference genes present in the CAZy database (Lombard et al., 2014). A first comparison was accomplished at strain level using the genome information from B. uniformis strains listed in Table 1. Consequently, we have recovered a total of 622, 626, 605, 624, and 671 CAZy genes for B. uniformis CECT 7771, B. uniformis ATCC 8492, B. uniformis CL03T00C23, B uniformis dnLKV2, and B. uniformis 3978-T3i, respectively. Those genes were grouped into non-redundat CAZy family domains, then families were split in case of multi-domain proteins and finally compared via Venn diagrams to disclose function uniqueness in every single strain (Figure 3A). We found a similar number of families in all strains and the assembled B. uniformis CECT 7771 genome encodes 115 different CAZy functional families whereas other strains encode up to 120. Interestingly, the five strains share in average 88% (103/117) of CAZy families and up to 90% (122/135) in the case of pairwise comparisons. When strain-specific functions were explored, we found that all strains showed unique features encoded in their respective chromosomes (Table S3). Although the glycosidehydrolase (GH) function was predominantly recovered in a set of unique genes in each strain, the presence of glycosyl transferase (GT) activities was only observed in B. uniformis CECT 7771 and 3978-T3i strains. The GT catalyzes the transference of sugar moieties from an activated donor to a specific substrate in an enzyme-dependent manner. In addition to the GT genes exclusively present in the genomes of strains of B. uniformis, the strains the B. uniformis CECT 7771 and 3978-T3i strains exhibited the largest number of GT enzyme encoded genes (30) when compared with counterparts (25 genes per strain on average). Therefore, the gain of GT genes would be indicative of the higher capability of these strains to synthesizing disaccharides, oligosaccharides, polysaccharides, and even glycolipids. Accordingly, B. uniformis CECT 7771 could represent an important and unexplored reservoir of GT functions for production of natural compounds with an ample range of biotechnology applications.

On the other hand, comparisons at species level were done using similar analysis as the previous one but with the genomic information of the annotated genes in the CAZy database and available for Bacteroides species such as: B. xylanisolvens XB1A, B. cellulosilyticus WH2, B. thetaiotaomicron VPI5482, B. vulgatus ATCC 8482 and B. fragilis NCTC 9343. Strikingly, the number of CAZy genes encoded in other Bacteroides species is lower than that detected in B. uniformis strains (Figure 3B). The above observation reveals a new degree of specialization of B. uniformis strains that outnumber the CAZy genes found in other Bacteroides species even in B. thetaiotaomicron, the archetype for characterization of PULs and B. cellulosilyticus previously described to encode an expanded glycolytic potential (McNulty et al., 2013). In the same line of thinking, almost the 23% (27/115) of B. uniformis CECT 7771 CAZy genes cannot be identified in other species. Moreover, the number of CAZy genes present in Bacteroides spp. genomes seems to be indicative of their role as commensals or potential pathogens. Thus, in species like B. vulgatus and B. fragilis, which are considered as pathobionts (Wexler, 2007) the number of these genes are abnormally lower than in other species; this is particularly notable for B. fragilis which has almost half of genes found in B. uniformis strains.

Finally, a searching for potential antibiotic resistance and virulence genes was completed by comparing the full set of ORFs predicted to be encoded in the B. uniformis CECT 7771 genome according to the proteins annotated in the ARDB, MvirDB, and VFDB databases. We detected presence of a β-lactamase protein encoded in the ORF BUNIF7771_3570 (99% sequence identity against the AAA66962 ARDB entry) and presence of the TetQ protein encoded in the ORF BUNIF7771_4507 (90% sequence identity against the AAS83507 ARDB entry), both proteins with homologs in other B. uniformis genomes analyzed in above sections. An additional KEGG based functional analysis of ORFs predicted to be encoded by the B. uniformis CECT 7771 genome has revealed the presence of a total of 24 genes corresponding to 11 different multidrug efflux pumps, one D-alanyl-D-alanine dipeptidase (vanX, K08641) associated with vancomycin resistance, and a couple of genes associated with cationic antimicrobial peptide resistance.

In the light of the wide repertoire of genes dedicated to polysaccharide degradation identified in B. uniformis CECT 7771 genome, we wanted to investigate those expressed in response to complex carbon sources that can be present in the intestinal tract (GIT) such as mucin, a heavily glycosylated secretion product of the intestine epithelial cells. To do so, the ability of B. uniformis CECT 7771 to grow in the presence of glucose, gum arabic, WBE, inulin, mucin, or pectin as primary carbon source was evaluated. Growth patterns differed (Figure 4A) essentially regarding the log phase and the yield at stationary phase. To quantify the ability to growth in the different carbon sources we also calculated the doubling time during the exponential growth phase (Figure 4B). In general, we observed that all carbon sources except for WBE increased the doubling time of the bacterium compared to glucose; this effect was especially significant for gum arabic (p < 0.022). Notably, WBE showed no fitness cost and boosted the B. uniformis CECT 7771 growth reducing the doubling time from 91 to 62 min (p < 0.057). The opposite effects of gum arabic and WBE on growth fitness could be explained by their composition. Whereas gum arabic consist of a complex mixture of branched polymers of galactose (39–42%), rhamnose (12–16%), arabinose (24–27%), and glucuronic acid (15–16%) (Ali et al., 2009), the WBE is principally composed of AXOS (>69%), which are a mixture of low molecular weight xylo- and arabinoxylo-oligosaccharides (Cargill, Antwerp, Belgium). The complex polysaccharides contained in gum arabic, inulin, and pectin seem to be fermented with more difficulty by B. uniformis CECT 7771 probably by the lacking of the proper de-branching enzymes useful in the first hydrolytic stages from poly- to oligosaccharides. Therefore, B. uniformis CECT 7771 is more specialized in the utilization of oligosaccharides than in very complex carbohydrates, which could also be related to the strain origin (isolated from stools of a breast-fed infant).

With the aim to complement the study above, we also did an exploratory genome-wide expression analysis in the presence of gum arabic, WBE, inulin, mucin, and pectin. This analysis was done using a RNAseq approach with RNA pools from respective replicates obtained in different and independent experiments. Using comparative expression algorithms implemented in GFOLD (Feng et al., 2012) we could discern a set of genes more probable to be differentially regulated under the respective conditions. By setting an increased threshold for selection (see Methods) we obtained a total 633 genes with a high probability to be differentially expressed among all conditions (Figure 5A). Interestingly, the gum arabic and mucin seem to induce and attenuate the expression of a larger set of genes in B. uniformis CECT 7771. Moreover, both carbon sources similarly modify the expression of 152 out of 633 (24%) genes differentially expressed in all conditions inducing a similar response. This is consistent with what graphically retrieved when a hierarchical clustering approach was used to study the particular gene expression patterns for the 633 genes with plausible differential expression (Figure 5B). Globally, the exposure of B. uniformis CECT 7771 to gum arabic and mucin appears to have a deep impact generating more transcriptional activity across its genome conversely to what observed in the response to WBE and inulin, which trend to trigger more concrete expression signals confined to a few and well-delimited gene clusters (Figure S1).

Despite the consistency among the metabolic functions up- and down-regulated depending on the carbon source used by B. uniformis CECT 7771 for growing, we validated the changes in expression of a set of 8 genes by qPCR, using replicated and non-pooled total RNA samples remaining from the exploratory whole-transcriptome experiment. These genes were selected for their possible role in the host physiology. In particular we analyzed the relative expression of following genes: BUNIF7771_0544 (K01580—glutamate decarboxylase), BUNIF7771_0548 (K20265—glutamate:GABA antiporter), BUNIF7771_1883 (K00929—butyrate kinase), BUNIF7771_0387 (K09789—pimeloyl-ACP methyl ester esterase), BUNIF7771_3732 (K01915—glutamine synthase), BUNIF7771_4131 (K00208—enoyl-ACP reductase I), BUNIF7771_1668 (K01813—L-rhamnose isomerase), and BUNIF7771_3473 (K01805—xylose isomerase). When we compared the respective and relative gene expression to that found in the presence of glucose, we confirmed almost all the gene expression patterns inferred from the transcriptome analysis (Figure S2).

Once we corroborated the expression pattern observed for B. uniformis CECT 7771 genes probably involved in GABA production by qPCR, we further analyzed the GABA concentration in supernatants of cultures supplemented with mucin and pectin that showed a clear over-expression of glutamate decarboxylase (BUNIF7771_0544) and glutamate:GABA antiporter (BUNIF7771_0544) genes (Figure S2). We observed that mucin and pectin increased the GABA concentration in the extracellular media by 45 and 63%, respectively, compared to glucose (Figure 6). These results together with the over-expression pattern inferred from the RNAseq approach, definitively confirm the role of different carbon sources in production of bioactive derived metabolites other than SCFAs, and, on the basis, suggest that innovative synbiotic products combining B. uniformis CECT 7771 with pectin and/or glycans similar to those covalently attached to mucins could be developed and functionally tested in pre-clinical trials.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.