Genome Structure of the Symbiont Bifidobacterium pseudocatenulatum CECT 7765 and Gene Expression Profiling in Response to Lactulose-Derived Oligosaccharides

Bifidobacterium pseudocatenulatum CECT 7765 was isolated from stools of a breast-fed infant. Although, this strain is generally considered an adult-type bifidobacterial species, it has also been shown to have pre-clinical efficacy in obesity models. In order to understand the molecular basis of its adaptation to complex carbohydrates and improve its potential functionality, we have analyzed its genome and transcriptome, as well as its metabolic output when growing in galacto-oligosaccharides derived from lactulose (GOS-Lu) as carbon source. B. pseudocatenulatum CECT 7765 shows strain-specific genome regions, including a great diversity of sugar metabolic-related genes. A preliminary and exploratory transcriptome analysis suggests candidate over-expression of several genes coding for sugar transporters and permeases; furthermore, five out of seven beta-galactosidases identified in the genome could be activated in response to GOS-Lu exposure. Here, we also propose that a specific gene cluster is involved in controlling the import and hydrolysis of certain di- and tri-saccharides, which seemed to be those primarily taken-up by the bifidobacterial strain. This was discerned from mass spectrometry-based quantification of different saccharide fractions of culture supernatants. Our results confirm that the expression of genes involved in sugar transport and metabolism and in the synthesis of leucine, an amino acid with a key role in glucose and energy homeostasis, was up-regulated by GOS-Lu. This was done using qPCR in addition to the exploratory information derived from the single-replicated RNAseq approach, together with the functional annotation of genes predicted to be encoded in the B. pseudocatenulatum CETC 7765 genome.


INTRODUCTION
Indigenous intestinal bacteria are known to be equipped with an array of genes coding for uptake systems and complex enzymatic machinery that facilitates the utilization of oligo-and polysaccharides. This constitutes a major adaptation mechanism to the main energy sources available in the large intestine for symbiotic bacteria (Jost et al., 2015). The role of breastfeeding in defining the composition of the gut microbiota, characterized by the dominance of bifidobacteria in infants, constitutes the best example of this adaptation process (Penders et al., 2006;Palma et al., 2012;Olivares et al., 2014). Human milk oligosaccharides are highly diverse and complex glycans, which seem to have evolved naturally, shaping the species that inhabit the infant gut (Koropatkin et al., 2012). Thus, species of the genus Bifidobacterium that predominate in breastfed babies (B. longum subsp. infantis and B. bifidum, the so-called infant-type bifidobacteria) are known to possess the genetic and protein machinery (oligosaccharide binding proteins, fucosidases, lacto-N biosidase, etc.) necessary to utilize humanmilk oligosaccharides. This confers a competitive advantage to bifidobacteria, whereby they outnumber other intestinal bacteria (Yoshida et al., 2012;Garrido et al., 2013;Viborg et al., 2014). This adaptation to diet may, in turn, also define the symbiotic host-microbe interactions and their biological role in human health (Sanz, 2015). Epidemiological studies have demonstrated that breast-feeding confers benefits for early and long-term health, improving intestinal transit and reducing the incidence of infections and non-communicable diseases (e.g., obesity and type-2 diabetes; reviewed by Sanz, 2015). It is not easy to confirm whether these effects are a direct consequence of the humanmilk-induced microbiota pattern; however, it is biologically plausible that bifidobacteria play a potential role, which is supported to some extent by existing clinical (Braegger et al., 2011;Miller and Ouwehand, 2013;Sanz, 2015) and pre-clinical data (Cano et al., 2013;Hayes et al., 2014;Moya-Perez et al., 2014, 2015Reichold et al., 2014;Elian et al., 2015;Srutkova et al., 2015). Consequently, alternative oligosaccharides have been produced to try to mimic the role of human milk oligosaccharides in the infant's microbiota, as well as to exert "bifidogenic" effects in adults, in the form of food ingredients or supplements (Marriage et al., 2015). Among these, long-chain inulin-type fructans (FOS) and short-chain galacto-oligosaccharides (GOS) are the most commonly used, particularly in infant formula attempting to provide the beneficial effects of breast-milk (Oozeer et al., 2013). Nevertheless, further research is required to develop oligosaccharides that more closely resemble their natural counterparts, and to shed light on their interactions with specific indigenous bacteria, as well as their health consequences.
Studies in vivo reveal that lactulose-derived GOS (GOS-Lu) have higher resistance to gastrointestinal digestion and less absorption in the small intestine of rats than GOS derived from lactose. These properties were attributed to the strong resistance of galactosyl-fructoses to the hydrolytic action of mammalian digestive enzymes (Hernandez-Hernandez et al., 2012b). Conventional GOS (derived from lactose) and GOS-Lu bear different structural features based not only on monosaccharide composition but also on the degree of polymerization, anomeric configuration, isomer composition and types of glycosidic linkage. Specifically, the predominant glycosidic linkage type in GOS-Lu is β-(1→6), whereas the main glycosidic linkage present in commercial GOS is β-(1→4) for instance Vivinal R (Friesland Campina, The Netherlands). Additionally, GOS-Lu also contains galactosyl-fructoses and galactobioses containing 1→2 and 1→5 linkages, which are not found in conventional GOS. These differences in glycosidic linkage types are thought to be crucial in the potential bioactivity of GOS-Lu type oligosaccharides.
In this study, we have investigated the ability of Bifidobacterium pseudocatenulatum CECT 7765, isolated from stools of a breast-fed infant, to utilize galacto-oligosaccharides. This is the first attempt to understand its origin, as this is generally considered to be an adult-type bifidobacterium, and to improve its functionality as a potential probiotic, bearing in mind this strain has proven pre-clinical efficacy in obesity and cirrhosis experimental models (Cano et al., 2013;Moya-Perez et al., 2014, 2015. In addition, we have evaluated the potential advantage of using a GOS-Lu instead of a GOS because of its persistence in the intestine and slower fermentation in proximal colon Martinez-Villaluenga et al., 2008;Hernandez-Hernandez et al., 2012b), as well as its capacity to selectively stimulate the growth and/or activity of bifidobacterial species (Marin-Manzano et al., 2013). Previous research also shows GOS-Lu are effective in a rat model of experimental colitis, presumably by exerting immunomodulatory effects associated with increased short-chain fatty-acid production (Algieri et al., 2014). Additionally, the presence of nontransgalactosylated lactulose, a prebiotic, instead of lactose in the GOS-Lu mixture could also provide additional benefits such as a lower calorific content than conventional GOS or other beneficial properties attributed to lactulose Martinez-Villaluenga et al., 2008).
To address this study, we have investigated the molecular response of B. pseudocatenulatum CECT 7765 when cultured in the presence of either glucose or GOS-Lu as carbon source. To do so, we have used three different high-throughput approaches: (i) genomic DNA sequencing for whole-genome assembly and functional annotation of the strain studied; (ii) a preliminary and exploratory single-replicated transcriptome approach of RNA pools to detect potential signals of over-expression across the B. pseudocatenulatum CETC 7765 genome during GOS-Lu fermentation, followed by validation of certain gene expression patterns by qPCR; and (iii) metabolite analysis of GOS-Lu species using gas chromatography coupled to mass spectrometry (GC-MS) to understand the upregulated metabolic pathways, the final output and potential biological consequences.

Bacterial Growth
Bifidobacterium pseudocatenulatum CECT 7765 was isolated from a breast-fed infant subject of a prospective observational study carried out in a cohort of 164 healthy full-term newborns, with a first degree relative affected by celiac disease. Ethics committee approvals for that study, according to the Helsinki Declaration of 1983, are already published (Palma et al., 2012). B. pseudocatenulatum CECT 7765 was grown overnight in low glucose modified MRS broth (10 g/L peptone, 8 g/L meat extract, 4 g/L yeast extract, 5 g/L sodium acetate, 2 g/L di-ammonium citrate, 0.2 g/L magnesium sulfate, 0.05 g/L manganesum sulfate, 2 g/L di-potassium phosphate, 0.1% v/v polysorbate 80) supplemented with 0.05% (w/v) L-cysteine and 0.5% (w/v) glucose at 37 • C under anaerobic conditions into Whitley DG250 Anaerobic Workstation (don Whitley Scientific, Inc., Shipley, UK). A 5 mL aliquot of an overnight culture was obtained for genomic DNA isolation. For RNA isolation, over-day cultures were obtained by refreshing 1/100 overnight cultures in pre-warmed and oxygen-depleted modified MRS media supplemented with 0.05% (w/v) L-cysteine and containing 1% (w/v) glucose or 1% (w/v) GOS-Lu as sole carbon sources. Early exponential growth phase cultures were collected after 7-8 h of incubation, when optical density at 600 nm (OD 600) was approximately 0.4, to optimally detect differential gene expression patterns, as previously described (Garrido et al., 2012), and to measure the residual fraction of GOS-Lu components. Stationary growth phase cultures were collected to measure accumulation of branched-chain amino acids in cell-free culture supernatants by LC with fluorescence detection. In both cases, cells were pelleted by centrifuging at 4 • C and 2,500 × g for 20 min, and supernatants were aspired off and filtered using 0.22 µm disposable filters (Millipore).

Nucleic Acid Isolation
DNA and RNA from respective bacterial cultures were isolated using MasterPure TM Gram Positive DNA Purification Kit (Epicentre) with slight variations over manufacturer's instructions. Briefly, a cell lysis step was improved by incubating cell suspension with 500 µg Lysozyme (Sigma, Cat #62970) and 20 U Mutanolysin (Sigma, Cat #M9901) for 60 min at 37 • C. For DNA isolation, samples were incubated with RNase A at 37 • C for 60 min, whereas for RNA isolation samples were incubated with 2 U DNase I (Epicentre) at 37 • C for 60 min instead of the RNase A treatment.

High-Throughput Sequencing
Ten µg genomic DNA from B. pseudocatenulatum CECT 7765 were sent to Eurofins Genomics GmbH (Ebersberg, Germany) to produce a shotgun library by fragmentation and end repair of DNA with an insert size of 300-400 bp and 2 bp × 150 bp (paired-end) configuration. The preliminary and exploratory transcriptome analysis was achieved by pooling RNA from three independent experiments in an equimolar mixture of glucose-or GOS-Lu-derived samples, respectively. Then, 30 µg total RNA per pool were also sent to Eurofins Genomics GmbH (Ebersberg, Germany) to produce cDNA libraries with an insert size of 150-400 bp and prior rRNA depletion using RiboZero TM Magnetic Kit Gram-Positive Bacteria (Epicentre). DNA and cDNA libraries were pooled and sequenced in one MiSeq cell flow allowing 1/10 proportion of DNA library against cDNA libraries.

Data Analysis
The B. pseudocatenulatum CECT 7765 genome was assembled using the MIRA assembler (Chevreux et al., 1999). Scaffolding of contigs was assisted by SSPACE (Boetzer et al., 2011) and scaffold reordering was predicted by using comparative genomics and whole-genome alignment algorithms implemented in MAUVE (Darling et al., 2010) and draft genomes of close species such as B. pseudocatenulatum DSM 20438 (Accession number NZ_ABXX02000001). Predicted joints were corroborated by PCR and Sanger sequencing. Gene prediction and functional annotation were performed using tRNAscan-SE (Lowe and Eddy, 1997), RNAmmer (Lagesen et al., 2007), Prodigal (Hyatt et al., 2010), KEGG Automatic Annotation System (Moriya et al., 2007), SMART database (Letunic et al., 2012), Pfam database (Finn et al., 2014), CAZy database (Lombard et al., 2014), CAT server (Park et al., 2010), ScanProsite (de Castro et al., 2006), and Artemis (Rutherford et al., 2000). Sequence information supporting the B. pseudocatenulatum CECT 7765 genome assembly was submitted to the European Nucleotide Archive (ENA) where it is publicly available under primary accession number PRJEB6926. Annotated 5S, 16S, and 23S rRNA gene sequences from B. pseudocatenulatum CECT 7765 are publicly available under ENA accession numbers LN624223, LN624224, and LN624525, respectively. Comparative genomics among Bifidobacterium species was accomplished using BRIG (Alikhan et al., 2011) and available and complete genome information from B. breve UCC2003 (NC_020517.1), B. dentium Bd1 (NC_013714.1), B. longum subsp. infantis ATCC 15697 (NC_017219.1), and B. pseudocatenulatum DSM 20438 (NZ_ABXX02000001, draft genome). The quality filtering and trimming of glucose-and GOS-Lu-derived RNA-seq analysis was performed using FASTX-toolkit 1 . Read mapping was assisted using the local alignment Blast algorithm (Altschul et al., 1990) and selecting alignments >50% of read length (>70 nt) and 100% identity. Read counts were normalized using RPKM (Mortazavi et al., 2008) and the exploratory differential expression among glucose-and GOS-Lu-derived transcriptomes was measured with GFOLD. This analysis tool is able to detect potential trends in gene expression in unreplicated data, requiring further evaluation by conclusive methods like qPCR (Feng et al., 2012). In order to increase the stringency for detecting plausible signals of differential expression, we only selected genes with GFOLD score ≤ −1 or ≥1. Sequence information supporting the B. pseudocatenulatum CECT 7765 transcriptome analysis was submitted to the ENA where it is publicly available under primary accession number PRJEB6928.

Quantitative PCR
The genes BPSEU7765_0088, BPSEU7765_0773, BPSEU7765_0523, BPSEU7765_0525, and BPSEU7765_1462 were selected from the preliminary and exploratory RNA-seq analysis to assess specific changes in expression by qPCR. The gene-specific oligonucleotides used for this aim are presented in the Supplementary 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 35 cycles of 10 s at 95 • , 20 s at 65 • , 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 GOS-Lu treatments.

Quantitative Analysis of GOS-Lu Consumption by GC-MS
Cell-free supernatants of B. pseudocatenulatum CECT 7765 cultures supplemented with GOS-Lu to replace glucose were obtained and analyzed by GC-MS using a two-step derivatization procedure (oximation and trimethylsilylation) according to previous methods . Trimethylsilyloximes (TMSO) derivatives of carbohydrates were identified by comparison of mass spectra and retention indices with standard derivatized carbohydrates, as described in a previous study (Hernandez-Hernandez et al., 2012a). Characteristic mass spectra and data previously reported in the literature were used to identify those carbohydrates unavailable as commercial standards. Carbohydrate quantitative data were obtained from GC-MS peak areas using the internal standard method. To do so, standard solutions from 0.003 to 1 mg of phenyl-β-D-glucoside, sucrose, and raffinose were prepared to calculate the corresponding response factors relative to internal standard and used to quantify mono-, di-, and tri-saccharides, respectively. Analytical standards of fructose, galactose, lactulose [β-D-galactopyranosyl-(1→4)-D-fructose], lactose

Quantitative Analysis of Branched-Chain Amino Acids
For valine (V), isoleucine (I) and leucine (L) analyses, samples were 1,000-fold diluted with 2 N acetic acid (BDH Prolabo) and subjected to an automatic pre-column derivatization with o-phthaldialdehyde (OPA; Sigma-Aldrich). For the derivatization step, 20 µl of diluted sample was mixed with 15 µl of 4 N NaOH (BDH Prolabo), 40 µl of OPA and 20 µl of 5% acetic acid in Milli-Q water and, then, 20 µl of the resulting mixture was injected into the LC system. Amino acids were separated by LC Gemini C-18 column (5 µm particle size, 250 mm × 4.6 mm i.d., Phenomenex) at a flow rate of 1 mL/min and 25 • C. Solvent A was 0.15 M anhydrous sodium acetate (Sigma-Aldrich)/HPLC grade methanol (BDH Prolabo;70:30,v:v) adjusted to pH 6.8 with glacial acetic acid (BDH Prolabo); solvent B was HPLC grade methanol/Milli-Q water (70:30, v:v). Elution was performed with a linear gradient as follows: 0-13 min, 50% B; 13-15 min, 100% B; 15-22 min, 100% B; 22-23 min, 50% B; 23-34 min, 50% B. Detection was performed by fluorescence using 340 and 455 nm for excitation and emission, respectively. Calibration curves of V, I, and L were built using commercial pure standards (0.1-1 ppm in 2 N acetic acid for I and L; 0.05-1 ppm in 2 N acetic acid for V) purchased from Sigma-Aldrich. Production of BCAAs was individually compared (V, I, or L) between treatments (GOS-Lu vs. Glucose) and differences were statistically analyzed with a one-sided t-test with Welch's correction.

B. pseudocatenulatum CECT 7765 Genome
The draft genome of B. pseudocatenulatum CECT 7765 comprised ∼2.25 Mbp containing six major super-scaffolds (ENA accession numbers: CDPW01000001 to CDPW01000006) with a 56.4% GC content, assembled from 138 contigs with N50 ∼145,000 bp and a theoretical coverage of 93X (MIRA assembler). At final assembly stage, the total number of genes predicted to be encoded by the B. pseudocatenulatum CECT 7765 genome was 1,879. Of these, 1,821 corresponded to coding genes, 54 tRNA genes, three rRNA genes clustered in a single operon, and one tRNA pseudogene with a CAT anticodon. This genome structure is quite similar to others from B. psudocatenulatum species recently sequenced (Alegria et al., 2014). Additionally, a Clustered Regularly Interspaced Short Palindromic Repeats (CRISPRs) region was predicted. This region is characterized by 57 repeats of the ATTTCAATCCACGCTCTCCATGAGGAGAGCGAC sequence and a CRISPR-associated gene (Cas) annotated as BPSEU7765_1382 gene immediately downstream of the  region. This feature indicates that B. pseudocatenulatum CECT 7765 can defend against bacteriophage attack with its innate immune system (Bhaya et al., 2011). Ribosomal RNA sequence information is publicly available in the ENA under accession numbers LN624223 to LN624225 for 5S, 16S, and 23S molecules, respectively. When the complete genomes (except for B. pseudocatenulatum DSM 20438) of representative species of the Bifidobacterium genus were compared, we could distinguish certain genomic regions differentially present in B. pseudocatenulatum CECT 7765 that are absent in other species, and even in one strain of the same species, B. pseudocatenulatum DSM 20438 (Figure 1). Consequently, we could distinguish six different genomic regions entirely and uniquely present in B. pseudocatenulatum CECT 7765, called Specificity Islands (SP1 to SP6) hereinafter. Functional analyses were made to disclose the potential gain-of-function in B. pseudocatenulatum CECT 7765 genome. In general terms, we found that functions distinctively found in the B. pseudocatenulatum CECT 7765 SPs were related to bacterial defense, and to carbohydrate transport and metabolism (Figure 1), and were sometimes duplicated. For instance, SP1 and SP2 show duplication in the cluster conformed by genes encoding DNA methylase, HTH_Tnp transposase and RelB antitoxin protein.
These duplication events can be explained by the potential presence of transposases, predicted to be encoded in the respective SPs.
As stated above, we observed a predominant gain of defense genes such as DNA methylases, restriction enzyme systems, and abortive phage infection proteins in the B. pseudocatenulatum CECT 7765 genome. Furthermore, there was a great variety of enzymes associated with carbohydrate metabolism and transport. Among them, we could distinguish genes coding for different metabolic functions such as glucosyltransferases, polysaccharide synthases, peptidoglycan-associated polymer synthases, and mono-(glucose and rhamnose) and tri-saccharide (raffinose) hydrolases, as well as multiple sugar transporters. Globally, the genome structure suggests that B. pseudocatenulatum CECT 7765 is strongly protected against phage infection by restriction modification systems located at SP2 (Figure 1), identified with locus tags BPSEU7765_0786 to BPSEU7765_0788, which account for a type I system with M, S, and R subunits, respectively. They appear to be additional to three other restriction modification systems, identified with locus tags BPSEU7765_708 to BPSEU7765_710 (type I), BPSEU7765_1106 (mrr protein of type IV system), and BPSEU7765_1110 to BPSEU7765_1112 (type III). Besides these bacterial immune system genes, we have also identified a CRISPR locus and a gene encoding a protein associated with the ability to abort active phage infections, the Abi protein (locus tag BPSEU7765_1115), previously reported in Lactococcus species (Anba et al., 1995;Garvey et al., 1995). Furthermore, B. pseudocatenulatum CECT 7765 seems to have a wide repertoire of genes involved in mono-, oligo-, and poly-saccharide metabolism, which could facilitate its survival in the large intestine where complex oligosaccharides are the main energy source. To better understand the potential of B. pseudocatenulatum CECT 7765 to metabolize a wide variety of oligo-and poly-saccharides, we have annotated all its genes encoding active enzymes of the carbohydrate metabolism, according to CAZy database (Lombard et al., 2014). Subsequently, we submitted the 1,821 ORF sequences, translated into amino acids, to the CAT server (Park et al., 2010). Thus, we have obtained 252 enzymes annotated with the CAZy system, encoded by the B. pseudocatenulatum CECT 7765 chromosome. A comparative analysis indicates that the CAZy enzymes present in B. pseudocatenulatum CECT 7765 outnumber those present in B. dentium Bd1 (128 CAZy enzymes), B. breve UCC2003 (83), B. longum ATCC 15697 (70), and B. pseudocatenulatum DSM 20438 (97), all of which are fully annotated in the CAZy database. We detected the specific group of enzymes exclusively present in B. pseudocatenulatum CECT 7765 by drawing Venn diagrams (Bardou et al., 2014) and disclosing the intersection lists of enzymes (Figure 2A). Consequently, we found that 33 CAZy families were exclusively found in B. pseudocatenulatum CECT 7765, whose associated functions are shown in Table 1. The above results indicate that B. pseudocatenulatum CECT 7765 should be able to grow in the presence of a wide variety of carbon sources. This assumption was reflected in the fact it thrived on media supplemented with GOS-Lu instead of glucose ( Figure 2B). Indeed, the doubling time (time to double the OD 600 absorbance) at the exponential growth phase was slightly lower in the presence of GOS-Lu than in the presence of glucose (73.3 ± 3.8 min vs. 76.1 ± 0.7, respectively, p = 0.48). The ease with which it utilizes this prebiotic substrate seems to be a particular feature of this strain, not exhibited by other Bifidobacterium strains when cultured with GOS as carbon source (Marcobal et al., 2010;Garrido et al., 2013;Watson et al., 2013). The species B. pseudocatenulatum has traditionally been considered an adulttype bifidobacteria, unlike B. breve and B. longum subsp. infantis, which are considered infant-type bifidobacteria (Pozo-Rubio et al., 2011). Nevertheless, B. pseudocatenulatum CECT 7765 was originally isolated from the stools of healthy breast-fed infants (Cano et al., 2013). The ability of this strain to occupy this niche and out-compete other colonizers could be explained by its adaptation to utilize a wide variety of oligosaccharides.
Although, breast-milk composition is much more complex and oligo-galactose as such has only been found in small amounts in human milk (Scientific Committee on Food, 2003), human-milk oligosaccharides could promote growth in vivo similarly to that observed in vitro when glucose is replaced by GOS-Lu.

Exploratory Analysis of the GOS-Lu-Associated Transcriptome
An exploratory RNA-seq approach has been used to analyze trends in differential gene expression in B. pseudocatenulatum CECT 7765 genome in response to either glucose or GOS-Lu as carbon source during anaerobic fermentation. With this approach, we detected transcripts in more than 99.3% of the coding genes initially predicted to be encoded by the B. pseudocatenulatum CECT 7765 genome, leaving sequence reads for just 12 genes unmapped. Globally, we wanted to explore potential signals of differential expression on the basis of measuring the GOS-Lu/glucose ratio of transcripts, thereby pinpointing the likely genomic regions responsible for GOS-Lu uptake and metabolism (Figure 3). This exploratory analysis showed that 76 different genes had a trend of upregulation (GFOLD score ≥ 1) when GOS-Lu was used Frontiers in Microbiology | www.frontiersin.org as carbon source ( Table 2) whereas 25 genes exhibited a tendency toward down-regulated (GFOLD score ≤ −1) under the same condition (Table 3). Altogether, we found five different gene clusters probably associated with GOS-Lu fermentation (see gene tags at top of Figure 3), which may represent a molecular signature for bacteria able to metabolize this carbon source. Interestingly, among these genes, there were membrane and periplasmic permeases as well as ABC transporters, beta-galactosidases, oligosaccharide metabolic enzymes, and transcriptional regulators, potentially associated with GOS-Lu fermentation (see functional annotation in Table 2). We found that B. pseudocatenulatum CECT 7765 genome encodes a total of seven different beta-galactosidases (BPSEU7765_0525, BPSEU7765_1410, BPSEU7765_1435, BPSEU7765_1462, BPSEU7765_1517, BPSEU7765_1518, and BPSEU7765_1737), and five of these showed and indication to be over-expressed in the presence of GOS-Lu. Conversely, the BPSEU7765_1410 gene showed a signal for down-regulation (−1.82) and BPSEU7765_1737 had a neutral score for potential up-or down-regulation. This would suggest this couple of encoded enzymes may not be specific for GOS-Lu utilization. In addition to these likely expression patterns for beta-galactosidase genes, at least 20 different transporters/permeases seemed to be involved in GOS-Lu uptake ( Table 2). Particularly, the BPSEU7765_0523 and BPSEU7765_0524 ABC-type multiple sugar permeases showed an important over-expression trend under GOS-Lu exposure as compared to glucose in the medium. The cluster of genes BPSEU7765_0522 to BPSEU7765_0527 showed the strongest indication of over-expression in the presence of the prebiotic tested (Figure 3; Table 2) and, therefore, they may potentially be the main elements mediating GOS-Lu uptake and utilization. This cluster would be primarily regulated by the BPSEU7765_0522 gene encoding for a LacI-type transcriptional regulator, which would simultaneously control expression of sugar transporters/importers (BPSEU7765_0523, BPSEU7765_0524, and BPSEU7765_0527) and the glycosyl hydrolases (BPSEU7765_0525 and BPSEU7765_0526).
As a result of the above preliminary and exploratory RNAseq analysis, five genes were selected to confirm over-expression. These were the BPSEU7765_0088 and BPSEU7765_0773 genes, which would participate in the BCAA metabolism according to KEGG annotation; and the BPSEU7765_0523, BPSEU7765_0525, and BPSEU7765_1462 genes, related to sugar metabolism (Tables 1 and 2). In all cases, their up-regulation was confirmed by qPCR experiments indicating a high expression level associated with GOS-Lu consumption (Figure 4). Linear regression analysis with GFOLD scores against RQ values for those genes tested indicated that there is a good correlation among results from both type of analyses (Pearson's r = 0.73).

Metabolic Processes Likely Activated by GOS-Lu
Based on the functional annotation of the genes potentially over-expressed with GOS-Lu as carbon source, according to the exploratory RNA-seq approach, we identified the metabolic pathways presumably activated during its fermentation using the KEGG module-based annotation. We observed fourteen different pathways/modules represented by the genes exhibiting a likely up-regulation. As expected, we detected probable activation of genes involved in galactose degradation (KEGG Module M00632) and saccharide, polyol, and in lipid transport (M00196, M00199, M00491, and M00207 modules). When GOS-Lu was used as carbon source we detected a tendency of overexpression in a wide variety of oxidoreductases, which would control and protect against oxidative stress by regulating NAD+, NADP+, or H 2 O 2 levels in B. pseudocatenulatum CECT 7765. Likewise, we also detected over-expression signals for several genes responsible for DNA repair, which probably counteract potential DNA damage from ROS (oxygen-reactive species) produced by central metabolism and metabolic byproducts. Additionally, we observed that GOS-Lu would promote the expression of genes coding for the biosynthesis of purine and pyrimidines (M00049, M00050, and M00053 modules) and branched-chain amino acid (BCAAs = V, L, and I; M00019 and M00570 modules). Moreover, we detected a probable overexpression of pyruvate kinase (K00873), a key component of the machinery for carbohydrate degradation which is responsible for phosphoenolpyruvate (PEP) production, the common precursor of reduction pathways that produce short-chain fatty acids (SCFAs) (den Besten et al., 2013). In order to confirm our preliminary findings, further supported by qPCR results showing a strong over-expression in the BPSEU7765_0088 and BPSEU7765_0773 genes associated to BCAA metabolism, we also measured the metabolic output of the BCAA metabolic pathways over-expressed during GOS-Lu fermentation. We measured BCAA released to the extracellular media during growth in the presence of either glucose or GOS-Lu, using LC and OPAderivatization (Table 4). We found an increase in all BCAA in the cell-free supernatants of B. pseudocatenulatum CECT 7765 cultures supplemented with GOS-Lu as carbon source, but only differences in leucine concentrations were significant (p ≤ 0.0371). Leucine concentrations were more than 11% higher in the supernatants of GOS-Lu cultures than in the supernatants of glucose cultures. B. pseudocatenulatum CECT 7765 was shown to ameliorate the metabolic and immune dysfunction of diet-induced obesity in mice, partly by reducing lipid absorption and exerting an anti-inflammatory effect (Cano et al., 2013;Moya-Perez et al., 2015). According to our present findings, the combination of this strain with GOS-Lu could theoretically mediate another beneficial effect against obesity via generation of leucine, given the role of this amino acid as nutrient sensor and inducer of satiety (Potier et al., 2009;Schwartz, 2013).
Finally, some genes presented in the SP5 also exhibited an over-expression sign as a result of GOS-Lu fermentation. Thus, the corresponding predicted ORFs BPSEU7765_1612 and BPSEU7765_1613 would appear to be associated with GOS-Lu consumption by B. pseudocatenulatum CECT 7765. Although no functional information could be inferred for those genes according to databases, their genomic context suggests their involvement in sugar-nucleotide metabolism, given that they are flanked by ABC-MS multiple sugar transporters, a UPD glucose  The GFOLD analysis can detect potential trends of gene expression in unreplicated data but require further evaluation by other methods such as qPCR. A GFOLD score other than zero indicates probable differential expression between two genes under two different conditions, an inference based on Bayesian probabilities of log 2 fold-change per gene across the genome (p = 0.01; Feng et al., 2012). Threshold for selection of probable up-regulated genes was ≥1 whereas GFOLD values ≤1 were set for selecting genes likely down-regulated during GOS-Lu fermentation.

GOS-Lu Species Preferentially Consumed by B. pseudocatenulatum CECT 7765
To further integrate metabolic data resulting from GOS-Lu fermentation by B. pseudocatenulatum CECT 7765, we measured the concentrations of derivated mono-, di-, and trisaccharides at baseline and after growth in the presence of GOS-Lu. We observed preferential utilization of saccharide species and, particularly, the disaccharide fraction. Specifically, consumption of disaccharides was estimated at 70% and the trisaccharide fraction decreased by about 37% on average (Table 5).
Additionally, we found that the extracellular concentration of monosaccharide components of lactulose increased by 30-40% (fructose and galactose, respectively; Table 5). This observation can be explained by presence of glycosyl hydrolases anchored to the cell-wall surface. To test this, we searched for an amino acid motif in the 252 CAZy enzymes detected in B. pseudocatenulatum CECT 7765 (Table 1) using the ScanProsite server (de Castro et al., 2006) and the amino acid patterns recognized by sortase enzymes, responsible for covalent attachment of proteins to cell-wall surface in Gram-positive bacteria. We did not find any enzymes harboring the amino acid pattern associated with sortase B activity N-P-[QK]-T-N, but we did find five enzymes containing the extended amino acid pattern recognized by sortase A proteins [LPSN]-[PAG]-X-T-G (Ton-That et al., 2004;van Leeuwen et al., 2014). Among these, only one harbored a glycosyl hydrolase domain, encoded by the gene BPSEU7765_0825 (S-A-I-T-G motif at C-ter). Although, the above analysis could explain the increased concentration of monosaccharides in the extracellular medium, we cannot discount the potential release of intracellular glycosylhydrolases to the culture medium due to spontaneous cell lysis during culture and/or supernatant preparation. The main GOS-Lu species largely consumed by B. pseudocatenulatum CECT 7765 were lactulose, 1,4-galactosyl-[1, 1-galactosyl]-fructose and 6 galactosyl-lactulose, which decreased by 71, 47, and 31%, respectively, after incubation. FIGURE 4 | Gene expression by qPCR. The data derived from preliminary and exploratory RNA-seq approach for genes BPSEU7765_0088, BPSEU7765_0523, BPSEU7765_0525, BPSEU7765_0773, and BPSEU7765_1462 was used as starting point to evaluate their over-expression by relative quantification in a qPCR assay. RQ ± SEM values are presented for all genes analyzed using three independent replicates per treatment and the average expression under glucose exposure as reference (see RQ ∼ 1 for these samples). In all cases the differential expression was significantly higher in samples from cultures with GOS-Lu (RQ > 5) than in those with glucose (p ≤ 0.05).

DISCUSSION
We have assembled the draft genome of B. pseudocatenulatum CECT 7765 using DNA sequence analysis of high-throughput sequencing data. This bifidobacterial strain was isolated from a breast-fed infant, and shows preclinical efficacy preventing obesity and metabolic dysfunction. Through comparative genomics we have detected certain strain-specific genome regions in B. pseudocatenulatum CECT 7765 indicating a gain-of-function associated with carbohydrate uptake and metabolism. Further features found in the B. pseudocatenulatum CECT 7765 chromosome are indicative of its genome stability, for example regarding protection against phage infections. These traits are generally considered relevant for potential probiotic applications. In particular, we have disclosed four restriction modification systems, a CRISPR system, and a system to abort phage infections. To our knowledge, bifidobacteria usually present up to three restriction modification systems (O'Connell Motherway et al., 2009, which means that B. pseudocatenulatum CECT 7765 would be the first bifidobacteria harboring a larger collection of such defense mechanisms. Regarding the number of enzymes related to saccharide metabolism, the CAZy annotation system (Lombard et al., 2014) showed us that B. pseudocatenulatum CECT 7765 contains a larger set of these proteins when compared with complete genomes of close species. Sixty-seven of the 252 CAZy enzymes detected in the B. pseudocatenulatum CECT 7765 genome were found to be exclusive to this strain when compared with the genetic information of close species, and they were grouped into 33 different CAZy families.  a Standard deviation in parenthesis (n = 3). b Reducing carbohydrates gives rise to two peaks corresponding to the TMS oximes of syn (E) and anti (Z) isomers after derivatization. c Reduction or increase in GOS-Lu mono-, di-, and tri-saccharide species was calculated as the relative percentage of the final concentration of those saccharides (GOS-Lu post) compared to their initial concentration (GOS-Lu pre), respectively.
An exploratory RNA-seq analysis using pooled samples enabled stool identify potential genes responsible for GOS-Lu fermentation, encoding a wide variety of sugar transporters and permeases. Also, the expression of five out of seven beta-galactosidases present in the B. pseudocatenulatum CECT 7765 genome was identified to be likely associated with GOS-Lu fermentation. The over-expression trend observed for some of those genes was further assessed through a qPCR approach. As a result, we confirmed the over-expression of the BPSEU7765_0088, BPSEU7765_0773, BPSEU7765_0523, BPSEU7765_0525, and BPSEU7765_1462 genes in response to GOS-Lu consumption. This fact, and the strong correlation between GFOLD scores and RQ values obtained for the respective genes, would confirm the reliability of the findings from the exploratory transcriptome analysis, as a result of the GOS-Lu consumption by B. pseudocatenulatum CECT 7765. Particularly, we found a specific gene cluster (BPSEU7765_0522 to BPSEU7765_0528 genes) that was strongly over-expressed when GOS-Lu was used as carbon source. Taking into account the pattern of GOS-Lu species consumption, we hypothesized this gene cluster could be directly involved in controlling the import and hydrolysis of di-and tri-saccharides shown to be preferentially taken-up by B. pseudocatenulatum CECT 7765. In addition, the functions of the probable over-expressed genes have been mapped to the main bacterial metabolic pathways using the KEGG hierarchical classification. We found that GOS-Lu fermentation would activate galactose, saccharide and polyol degradation, lipid transport, antioxidant response, DNA repair processes, and purine/pyrimidine and BCAA biosynthesis. We corroborated the specific response of these genes to GOS-Lu, demonstrating that GOS-Lu fermentation boosts leucine synthesis and release, by directly analyzing the metabolic products generated and performing qPCR of some of the transcripts tentatively over-expressed in the exploratory RNAseq analysis. Therefore, the use of GOS-Lu could contribute to promoting B. pseudocatenulatum CECT 7765 growth in the gut and, additionally, the molecular and metabolic outputs obtained from these synbiotic interactions indicate its likely beneficial anti-obesity effects. These effects could be related to the role of leucine as nutrient sensor and inducer of satiety, via activation of the mTOR-S6K signaling pathway on the hypothalamus, which would promote expression of anorexigenic peptides (Schwartz, 2013). Nevertheless, the extent to which amino acids produced by intestinal bacteria can confer health benefits to the human host in vivo has yet to be demonstrated.