In Silico Probiogenomic Characterization of Lactobacillus delbrueckii subsp. lactis A4 Strain Isolated from an Armenian Honeybee Gut

Simple Summary Although L. delbrueckii strains have been isolated from different bee products and beehive niches to our knowledge, this is the first report of the isolation of the L. delbrueckii strain from the gut of a honeybee. Furthermore, the genomic analysis of the strain revealed two major aspects of its genomic constitution: (1) it possesses many genetic determinants attributing probiotic properties for the honeybees, and (2) it shows a reduction in size which is a typical adaptation of bacteria undergoing transformation to an endosymbiont. These two observations motivated us to hypothesize that in the case of L. delbrueckii subsp. lactis A4, we witness the first case of an L. delbrueckii strain evolution to a honeybee endosymbiont. Abstract A Lactobacillus delbrueckii ssp. lactis strain named A4, isolated from the gut of an Armenian honeybee, was subjected to a probiogenomic characterization because of its unusual origin. A whole-genome sequencing was performed, and the bioinformatic analysis of its genome revealed a reduction in the genome size and the number of the genes—a process typical for the adaptation to endosymbiotic conditions. Further analysis of the genome revealed that Lactobacillus delbrueckii ssp. lactis strain named A4 could play the role of a probiotic endosymbiont because of the presence of intact genetic sequences determining antioxidant properties, exopolysaccharides synthesis, adhesion properties, and biofilm formation, as well as an antagonistic activity against some pathogens which is not due to pH or bacteriocins production. Additionally, the genomic analysis revealed significant potential for stress tolerance, such as extreme pH, osmotic stress, and high temperature. To our knowledge, this is the first report of a potentially endosymbiotic Lactobacillus delbrueckii ssp. lactis strain adapted to and playing beneficial roles for its host.


Introduction
It is not easy to estimate precisely the economic value of honeybees as pollinators of crops and other plants. However, in any case, worldwide, it ranges up to dozens of billions of USD annually [1]. One of the main factors affecting the honeybee's welfare is the composition of their gastrointestinal microbiota, which is well characterized and comprises several phylotypes [2]. A crucial role in the honeybee's welfare is played by several core species and genera such as Snodgrassella alvi, Gilliamella apicola, Frischella perrara, and Commensalibacter sp. but also members of the genus Bifidobacterium and the former genus Lactobacillus [3]. Their representatives play crucial roles in pollen degradation and

Sampling and Maintenance of the Strain
The sampling was carried out in an apiary situated in the Ararat province in the town of Artashat (40 • 24 18"N, 44 • 34 35"E). The strain A4 was isolated from the gut of a young worker honeybee sampled on 25.08.2017 by enrichment of 10% skim milk with the particles of the honeybee's gut. First, 50 µL from the inoculated milk was mixed with 1 mL of sterile peptone water (PW). Next, 10-fold serial dilutions in PW were made up to 10-4, then 100 µL of the 10 −2 and 10 −3 dilutions were plated on Petri dishes with De Man, Rogosa, and Sharpe (MRS) 1.5% agar (Merck, Darmstadt, Germany). After incubation for 36-48 h at 30 • C, single colonies were picked up and inoculated in 3 mL MRS liquid broth and incubated at 30 • C for 24 h. After that period, ten isolates were plated on MRS 1.5% agar with a sterile loop by the agar streaking method and incubated at 30 • C for 24 h. This procedure was repeated thrice, and the cultures were microscopically monitored for contamination. Finally, stocks of skim milk were prepared from each isolate [15] and were kept at −70 • C.

DNA Isolation
Total DNA was isolated from 3mL 24-h liquid culture in MRS broth obtained by single colony inoculation. The DNA was extracted using the "Gram Plus & Yeast Genomic DNA Purification Kit" (EURx, Gdansk, Poland) according to the manufacturer's instructions. The final elution was performed in 70 µL of the kit's elution buffer. The quality of the isolated DNA was monitored on a 0.8% agarose gel in a TBE buffer system, while the concentration was determined on a Quantus™ fluorimeter (Promega, Madison, WI 53711 USA). The isolated DNA was further stored at −70 • C.

DNA Techniques
First, for preliminary species determination, total DNA isolated from the strain was sent to Macrogen (Seoul, Republic of Korea) for 16S ribosomal RNA gene sequencing by the Sanger chain termination method with the primers pair 27F/1492R [16] on the ABI 3730xl System. 50 µL of the isolated DNA was shipped on dry ice to BGI (Tai Po, N.T. Hong Kong) for whole-genome sequencing on the DNBseq platform [17] for DNB-SEQ PE100/PE150 sequencing with a generation of 1 Gb of data.

Bioinformatic Analyses
The forward and the reverse reads of the 16S rRNA gene sequencing were assembled with SeqMan™ II v.5.01 software (DNASTAR Inc., Madison, WI 53705, USA), and the assembled contig was uploaded to the GenBank of the NCBI. The results were analyzed online at the NCBI by the Nucleotide BLAST program (https://blast.ncbi.nlm.nih.gov/Blast.cgi, accessed on 1 February 2023) [18]. The quality of the raw data obtained from the DBBseq platform was monitored by the online tool FastQC v. 0.11.9 [19]. The contigs and scaffolds were assembled by the online assembler SPAdes v. 3.15.4 [20]. The quality of the assembled contigs was assessed by the online tool Quast v. 5.2.0 [21]. The contigs shorter than 200 bp were filtered by the Filter FASTA v. 1.9.1.0 [22] online software. The last four steps were performed on the Galaxy server (usegalaxy.eu, accessed on 3 February 2023). To confirm the species attribution, the average nucleotide identity of the assembled genome was evaluated with the ANI calculator [23] located on the EzBioCloud (ezbiocloud.net/tools/ani, accessed on 11 February 2023). The assembled contigs were uploaded to the GenBank of the NCBI. The annotation was performed by the NCBI Prokaryotic Genome Annotation Pipeline (PGAP) (https://www.ncbi.nlm.nih.gov/genome/annotation_prok/, accessed on 11 February 2023) [24][25][26]. The initial genotypic characterization was performed on the Center for Genomic Epidemiology (genomicepidemiology.org/services/, accessed on 11 February 2023) using the following tools: PlasmidFinder 2.1 (database version: 18 January 2023) [27][28][29] for the presence of plasmid replicons; ResFinder 4.1 (database version: 24 May 2022 and ResFinderFG 2.0 (database version: 30 June 2022 [27,[29][30][31] for the presence of antibiotics resistance genetic determinants; PathogenFinder 1.1 [27,32] for the presence of pathogenicity genetic determinants; and MobileElementFinder v. 1.0.3 (database version: 9 June 2020) [33] for identification of mobile genetic elements which could be related to antimicrobial resistance genes and virulence factors. The presence of bacteriocins' genetic determinants was checked with the online tool BAGEL4 (http://bagel4.molgenrug.nl/, accessed on 11 February 2023) [34]. The presence of the genetic determinants determining potential probiotic properties was performed using the NCBI's Sequence Set Browser (https://www.ncbi.nlm.nih.gov/Traces/wgs/, accessed on 15 February 2023), followed by BLAST analyses against the NCBI's database.

Genome Assembly, Annotation, and Characterization
The results from sequencing the 16S rRNA gene of isolate A4 are available at the NCBI's GenBank under the accession number OP784418.1. From the raw data of the next-generation whole-genome shotgun sequencing, 145 contigs were assembled. The assembly statistics are summarized in Table 1. The assembled contigs of the draft genome are available at the NCBI's GenBank under accession number JAQSVE000000000. Both BLAST of the 16S rRNA gene sequence and ANI of the assembled genome attributed the isolate to Lactobacillus delbrueckii ssp. lactis based on the 16S rRNA gene homology and the average nucleotide identity, respectively, compared with the NCBI databases. The results from the annotation are presented in Table 2. No plasmid origins of replication were detected by the PlasmidFinder tool. Neither were found genetic determinants for antibiotic resistance and factors of pathogenicity by the ResFinder, ResFinderFG, and PathogenFinder tools. The MobileElementFinder tool revealed the presence of an intact ISL6 insertion sequence belonging to the IS3 family (IS150 group) and seventeen 3 -or 5 -truncated mobile elements belonging to the ISLre2, ISL3, IS4, IS30, IS110, and IS256 families.

Genotypic Characterization of L. delbrueckii A4
The results of the BLAST analysis of the genes with the potential to determine probiotic properties are presented in Table 3. All the listed database hits have the same length with the query sequence, a query cover of 100%, and a percentage of identity of 100.00%.

Discussion
The analysis of the draft genome assembly performed with the Quast tool (Table 1) revealed that the assembly quality has the potential to guarantee correct analysis of L. delbrueckii ssp. lactis A4 genome, mainly because of the number of the contigs greater than 1000 bp, the N50, N90, L50, and L90 values. The calculated genome size based on all assembled contigs was 1,881,475 bp. This value was not significantly different from the values of the genome sizes calculated based on the total length and the contigs larger than 1000 bp, which once again confirmed the good quality of the assembly. The calculated GC content was 49.84%-a typical value for the species [35]. Surprisingly, when it was compared to other L. delbrueckii ssp. lactis strains, the A4 stood out with its clearly reduced genome size (Table 4). More interesting results were obtained from the annotation of the genome, which revealed a reduction not only in the genome size but also in the number of genes. L. delbrueckii ssp. lactis A4 possesses the least total number of genes when compared with the other 14 randomly chosen strains of the same species, the least number of pseudogenes, and, with one exception, the least number of RNA coding genes. This tendency is also partially visible when comparing the numbers of the protein-coding genes; still, four other strains possess a lesser number.
It is quite difficult to speculate why L. delbrueckii ssp. lactis A4 possesses a reduced genome. However, considering its unusual origin (to our knowledge, L. delbrueckii was never reported to be isolated from a honeybee gut), the hypothesis of some specialization to the honeybee gut ecological niche could not be excluded. It is observed and documented that the symbiosis between insects and bacteria leads to mutual adaptations. One aspect of these symbiotic interactions is the reduction of the genome size of the symbiotic bacteria [36,37]. This phenomenon also occurs in honeybees, as documented for Apibacter sp. [38], and the obligate endosymbionts Gillamella apicola and Snodgrassella alvi [39]. The distinctiveness of the L. delbrueckii ssp. lactis A4, as an example of a strain undergoing such adaptation, is further supported by the fact that, till now, to our knowledge, this is the only representative of the species isolated from a honeybee gut microbiota. We believe that this hypothesis is also supported by the fact that the reduction in the number of the genes affects mainly the RNA coding genes and pseudogenes, but to a minor extent, the number of the protein-coding genes, and thus keeping its "metabolic potential" which in the case of the species L. delbrueckii which lacks virulence factors and pathogenicity determinants, is beneficial for the host.
Even though lactobacilli (in the broader meaning of the former big genus before its reclassification in 2020 [8]) are considered a part of the honeybees' core endosymbionts [3], till now, L. delbrueckii strains have not been reported to be isolated from the honeybees. So, logically, if the hypothesis of the genome adaption is true, what would be the role of L. delbrueckii ssp. lactis A4 for the honeybees, so is it selectively maintained long enough for genomic changes to occur? A logical answer to this question is that it possesses some crucial probiotic properties for the host. So, the next step of this research was to make a probiogenomic characterization, a notion proposed by de Jesus et al. to describe the analysis of the genome for the presence of genetic determinants giving probiotic properties [14], which in honeybees could be mainly antioxidant properties, exopolysaccharides synthesis, stress tolerance, adhesion to the gut epithelium and biofilm formation, as well as antagonistic activity towards pathogens.
It is documented that in human and animal models, lactic acid bacteria (LAB) and particularly L. delbrueckii, exert a beneficial protective effect against reactive oxygen species (ROS) on the intestinal mucosa and epithelium [40]. The studies of some of the mechanisms determining these antioxidant activities in LAB and lactobacilli revealed that they involve the synthesis of several enzymes and proteins [41,42], as well as the synthesis of exopolysaccharides (EPS) [43]. Nineteen enzymes and proteins involved in the common oxidative stress resistance are listed in different Lactobacillus species [41]. Genetic determinants for 7 of them were found in L. delbrueckii ssp. lactis A4: NAD(P)H-dependent oxidoreductase, DNA starvation/stationary phase protection protein, DsbA family protein, NAD(P)/FADdependent oxidoreductase, thioredoxin family protein, thioredoxin-disulfide reductase and thioredoxin (Table 3). Oxidative resistance genes are rather scarce among lactobacilli [41], so finding such a number in only one strain should not be considered a coincidence.
Additional antioxidant activity can be attributed to the chelator activity of the EPS, which can sequester heavy metal ions [43]; however, they can exert additional beneficial properties on their own. Different gene products participate in exopolysaccharide synthesis [44][45][46]. Genetic determinants for 9 of them were found in L. delbrueckii ssp. lactis A4 genome: exopolysaccharide biosynthesis protein, NADP-dependent phosphogluconate dehydrogenase, phosphoglucosamine mutase, beta-phosphoglucomutase, UDP-glucosehexose-1-phosphate uridylyltransferase, UTP-glucose-1-phosphate uridylyltransferase GalU, bifunctional UDP-N-acetylglucosamine, diphosphorylase/glucosamine-1-phosphate N-acetyltransferase GlmU, oligosaccharide flippase family protein and flippase (Table 3). They strongly suggest an ability of EPS synthesis, which in turn can contribute to the already mentioned protection from heavy metals, but also by interfering with the adhesion of pathogens, as well as by exerting immunomodulatory properties and helping the biofilm formation [46].
Despite being related to EPS production, adhesion, and biofilm formation, they depend on different genetic determinants [47][48][49][50]. Eight genes encoding SLAP domain-containing proteins, D-alanyl-lipoteichoic acid biosynthesis protein DltB, D-alanyl-lipoteichoic acid biosynthesis protein DltD, elongation factor Tu, type I glyceraldehyde-3-phosphate dehydrogenase and triose-phosphate isomerase were found within the L. delbrueckii ssp. lactis A4 genome (Table 3). They could determine good intestinal adhesion properties, which in turn can result in competitiveness with pathogens [48] but also in immunomodulating and anti-inflammatory effects, as shown in in vitro studies [47]. Additionally, the strain probably possesses the ability for biofilm formation because three genes determine the synthesis of AI-2E family transporters and cell wall metabolism sensor histidine kinase WalK, involved in signaling and quorum sensing mechanisms [49,[51][52][53].
The inhibition of pathogens' growth could also be considered a probiotic property. In general, this property is attributed to the synthesis of bacteriocins or low molecular weight acids, which decrease the pH to levels intolerable to most pathogenic species. However, this is obviously not the case for the A4 strain because of the lack of bacteriocin genetic determi-nants and because it expresses this ability in neutralized conditions. Two alternatives are the production of H 2 O 2 and biogenic amines. Lactate oxidase and pyruvate oxidase are two enzymes whose activity can lead to the production of H 2 O 2 [54,55], and their genetic determinants were found within the A4 genome. Furthermore, genes encoding eight enzymes participating in the amino acids metabolism, which could also be involved in biogenic amines synthesis [56], were found: pyridoxal-dependent decarboxylase, orotidine-5'-phosphate decarboxylase, carboxymuconolactone decarboxylase family protein, diphosphomevalonate decarboxylase, diaminopimelate decarboxylase, putative ornithine decarboxylase, and bifunctional phosphopantothenoylcysteine decarboxylase/phosphopantothenate-cysteine ligase CoaBC (Table 3). However, for the moment, we do not have enough data on which of those two mechanisms is involved in the inhibitory activity, if it is due to both of them or some other unknown mechanism.
Finally, a probiotic strain to survive and be selectively maintained within the gut microbiome, especially in the case of the honeybees, should be stress tolerant and resistant to acidic and alkaline pH, osmotic stress, and heat. Extreme pH values are found in the different parts of the honeybee's gut, while the feeding could achieve osmotic stress with honey. In addition, the body temperature could sometimes rise dramatically during a flight on a warm sunny day. Many genes are involved in coping with these extrema [14], and some were found within the L. delbrueckii ssp. lactis A4 genome-those encoding orotidine-5'-phosphate decarboxylase, S-ribosylhomocysteine lyase, phosphopyruvate hydratase, several AAA family ATPases, Na + /H+ antiporter NhaC, two-component system regulatory protein YycI, nucleotide exchange factor GrpE, molecular chaperones DnaK and DnaJ, heat-inducible transcriptional repressor HrcA, chaperonin GroEL, oligoendopeptidase F, CTP synthase, glucosamine-6-phosphate deaminase, aquaporin family proteins, aldo/keto reductase, ATP-dependent Clp protease ATP-binding subunit ClpX and F0F1 ATP synthase subunit alpha, beta, gamma, epsilon, B and C (Table 3). These findings suggest strong stress tolerance against the listed factors.

Conclusions
The whole-genome sequencing of L. delbrueckii ssp. lactis A4 genome suggests that this strain is an unusual one, probably undergoing adaptation to the honeybee gut conditions and thus becoming part of the symbiont microflora-a process never documented before for this species. The probiogenomic bioinformatic analysis supports this hypothesis because it revealed genetic determinants for a different probiotic for the host properties. The strain is stress-tolerant-a prerequisite to survive and to be selectively maintained within the insect's microbiota.

Conflicts of Interest:
The authors declare no conflict of interest.