All publicly available genome sequences of B. animalis subsp. lactis (50 strains including the type strain DSM 10140) and B. animalis subsp. animalis (8 strains including the type strain ATCC 25527) were downloaded from the NCBI microbe genome database on the 21st of November 2019 (Sayers et al., 2019).

Subspecies identification was either obtained from previously published articles (Lugli et al., 2019) or performed by employing the rpoA and 16S ribosomal DNA sequence. A >98% identity to the type strain genes was used as threshold and the genes should furthermore be different from the type strain of a related subspecies, in this case B. animalis subsp. animalis, as shown through a phylogenetic tree (data not shown).

The sequence quality was assessed and sequences with an average coverage of ≥30× and a contig number below 120 were considered acceptable for phylogenetic analysis. The quality of the genomes was also evaluated by checking that the length of the sequenced genome corresponds with the expected length of the genome, based on the type strain (Milani et al., 2014).

Other bifidobacterial species, which have been shown to harbor tet(W) (Ammor et al., 2008; Wang et al., 2017) were also downloaded from the NCBI microbe genome database on the 21st of November 2019 and included B. longum (14 strains, type strain NCTC11818), B. thermophilum (6 strains, type strain DSM 20212), B. bifidum (11 strains, type strain ATCC 29521), B. pseudolongum (4 strains, type strain DSM 20099), B. pseudocatenulatum (3 strains, type strain DSM 20438) and B. breve (41 strains, type strain NCTC 11815). All tet(W) sequences from other genera where the gene have been described (Scott et al., 1997; Chopra and Roberts, 2001; Flórez et al., 2006; Kazimierczak et al., 2006; Ammor et al., 2008; Palmieri et al., 2011; Schröder et al., 2012) and shared identity to the tet(W) gene found in B. animalis subsp. lactis were also downloaded from NCBI on the 21st of November 2019.

ResFinder (Zankari et al., 2012), with a 80% identity threshold, was used to search for the presence of tet(W) in the examined genomes and the Rapid Annotation using Subsystems Technology (RAST) server with default settings was used to annotate the genomes. The annotated genomes were downloaded in GenBank format from the RAST server (Aziz et al., 2008; Overbeek et al., 2014) and imported to CLC Genomics Workbench 20 (Qiagen Bioinformatics, Aarhus, Denmark), where the presence of tet(W), its flanking genes and presence of mobile genetic elements was examined. tet(W) nucleotide and protein sequences was extracted from the annotated genomes for further phylogenetic analysis. GC content of tet(W) and other genes was assessed by employing the DNA/RNA GC Content Calculator at ENDMEMO (Endmemo, 2020).

The blastN tool available at ISFinder (Siguier et al., 2006) with default settings was used to determine the identity of the mobile genetic protein next to tet(W) in B. animalis subsp. lactis and its sequence was used to search for its presence in other genomic regions in the B. animalis subsp. lactis genomes, which was performed in CLC Genomics Workbench 20 (Qiagen Bioinformatics, Aarhus, Denmark).

The phylogenetic analysis of tet(W) included both the nucleotide and protein sequences from B. animalis subsp. lactis (Supplementary Table 1) and tet(W) genes found in other bifidobacterial species and other genera where the presence of tet(W) previously have been published (Table 1) (Scott et al., 1997; Chopra and Roberts, 2001; Flórez et al., 2006; Kazimierczak et al., 2006; Ammor et al., 2008; Palmieri et al., 2011; Schröder et al., 2012). The nucleotide and protein tet(W) sequences was either extracted from the annotated genomes or from NCBI (Sayers et al., 2019).

ClustalX2 (Larkin et al., 2007) was used to perform a pairwise multiple alignment of the tet(W) sequences (Higgins and Sharp, 1988) and BioEdit (Hall, 1999) was used to remove gaps and unpaired ends. The nucleotide phylogeny was built by evolutionary analysis by the Maximum Likelihood method and Tamura-Nei model by MEGA X (Tamura and Nei, 1993; Kumar et al., 2018) and the amino acid phylogeny was built by evolutionary analysis by Maximum Likelihood method and JTT matrix-based model also by MEGA X (Jones et al., 1992; Kumar et al., 2018). Number of single nucleotide polymorphisms (SNPs) and single amino acid polymorphisms (SAPs) was obtained from the multiple alignment output from MEGA X that was used to build the phylogenetic relationships.

The genomes, either fully assembled or contigs were annotated by Prokka, which annotates genomes through the use of different tools including Prodigal (coding sequences), RNAmmer (Ribosomal RNA genes), Aragorn (Transfer RNA genes), SignalP (Signal leader peptides) and Infernal (Non-coding RNA) (Seemann, 2014). Prokka annotation is a requirement for using Roary, since the .gff file (file containing sequences and annotations) provided by Prokka is used by Roary to create a multi-FASTA alignment of all the core genes (Page et al., 2015). Roary was set to perform nucleotide alignment using MAFFT and a Blastp percentage identity at 80% (Katoh, 2002). FastTree was used to produce an approximately maximum-likelihood phylogenetic tree from the core gene alignment file, which was visualized by MEGA X (Price et al., 2009, 2010; Kumar et al., 2018).

A total of 50 publicly available B. animalis subsp. lactis strains including the type strain DSM 10140 were downloaded from NCBI and consisted either of contigs or assembled genomes (Supplementary Table 1). The sequence quality was assessed and sequences with an average coverage of ≥30 fold and a contig number below 120 were considered acceptable. On this basis, six strains (B420, DS1_2, BI-04, IDCC4301, CF3_2, AD011) were excluded from the study. The genomes of CNCM I-2994 (Chervaux et al., 2011) and AD011 (Kim et al., 2009) had both been sequenced by Sanger shotgun sequencing and consist of complete genomes. However, AD011 has previously been shown to exhibit a poor sequence quality and was therefore excluded (Garrigues et al., 2010), CNCM I-2994 was not excluded from the study. A total of 44 genome sequences were therefore acceptable for further phylogenetic analysis.

The B. animalis subsp. lactis genomes exhibited a size of 1.91–2.08 Mb with a GC content of 60.0–60.6% (Supplementary Table 1), which is in agreement with data for the type strain of the subspecies (Milani et al., 2014).

Subspecies identification was either obtained from previously published articles (Lugli et al., 2019) or performed by analysis of the rpoA and 16S ribosomal DNA sequence.

The majority of the B. animalis subsp. lactis strains originated from human feces, but also from food samples, dietary supplements and domestic pigs, chimpanzees, rabbits, vervet monkeys, a barbary macaque, three different dog breeds and one strain, the genomic unique ATCC 27673 (Loquasto et al., 2013) originated from sewage (Supplementary Table 1). Species within the bifidobacterial genera are commonly found in the gastrointestinal tract of various animals, the human oral cavity and sewage (Milani et al., 2014) and the strains in this study therefore represent the most common habitats of bifidobacteria.

Since B. animalis subsp. lactis is included in a wide range of probiotics, it cannot be excluded that the strains isolated from human feces, domestic pigs and dogs originate from ingested products such as probiotics. However, the strain collection also include strains such as Bl12 that has been isolated from a healthy patient, which has not ingested probiotic products (Milani et al., 2013) and rabbits and monkeys have with high likelihood not been exposed to probiotics and these strains are therefore expected to be diverse from the industrially exploited strains. The genome sizes of the different strains also vary, which also indicate that the strains are diverse (Supplementary Table 1). Most of the strains are isolated or submitted to NCBI between year 2006–2018, which reflect the increased focus on probiotics in the last decades (Gogineni, 2013), while the type strain DSM 10140 originates from 1997 (Supplementary Table 1). However, the submission date of the genome sequences to NCBI does not necessarily reflect the time of isolation as some strains are isolated even earlier.

B. animalis subsp. lactis has previously been shown to be a strict monophyletic bifidobacterial taxon that has recently evolved (Milani et al., 2013), however, some diversity is observed between the strains within the subspecies based on the presence of truly unique genes in some of the strains (Lugli et al., 2019). The strains with the highest number of truly unique genes are also included in this study. It is therefore concluded that the strains included in the current study represent the diversity within the subspecies.

A 1920 bp tet(W) gene flanked by genes annotated as mobile element protein (966 bp), with inverted repeats at both ends of 50 bp and a hypothetical protein (HP) of unknown function (183 bp) was found in the majority of the studied B. animalis subsp. lactis strains (38 out of 44). These genes exhibit similar GC content (51.01–53.23%), which is lower than the flanking genes in the genetic region (52.46–62.25%) (Figure 1) and the average of the genome (60.0–60.6%) (Supplementary Table 1). tet(W) genes found in non-bifidobacterial and bifidobacterial species exhibit a GC content of 52.19–53.18%, indicating that tet(W) genes generally exhibit a GC content around 53%.

The three strains originating from dogs (2007B, 2010B, 2011B) did not encode tet(W), the mobile element protein or the HP (Figure 1 and Supplementary Table 1). Two strains (DS28_2, LMG P-17502_2) only encoded the tet(W) gene, while LMG P-17502 encoded tet(W) and the mobile element protein (Figure 1). UBBLa 70 exhibited a large deletion in the tet(W), with only 117 bp remaining and two strains (ATCC 27673, 1528B) encoded a truncated version of the mobile element protein. This indicate that the three genes have been present originally in B. animalis subsp. lactis but have been subject to deletion in some strains. Despite these differences, the presence of tet(W), the putative mobile element protein and the HP are highly conserved within B. animalis subsp. lactis strains. This conservation was even observed in the strains that are more genomic unique which include ATCC 27673 and 1528B, and the Bl12 strain and the strains isolated from monkeys and rabbits. This suggest that the genetic organization surrounding tet(W) is not only present in the industrially exploited B. animalis subsp. lactis strains.

The tet(W), the mobile element protein and the HP genes were positioned in the same genomic context in the majority of the examined strains, however, in a few strains, alterations downstream (DS28_2, LMG P-17502_1, LMG P-17502_2, 2007B, 2010B, 2011B) and upstream (2011B) (Figure 1) of the three genes were observed. These were the same strains that exhibited complete or partial deletions of the tet(W), the mobile element protein and HP genes.

The genomic position of tet(W) was also reported by Rozman et al. (2020). They suggest that tet(W) and its flanking genes from the HP before the IS element to the HP after isochorismate pyruvate-lyase (Figure 1), based on nucleotide bias and codon usage bias, is part of a putative genomic island that has co-evolved together with B. animalis subsp. lactis and originate from an ancestral host (Guo et al., 2012; Bertelli et al., 2017). The codon usage bias corresponds with the gene GC content being lower in these genes compared to the rest of the genome. Genomic islands are defined as clusters of genes in bacterial genomes of probable horizontal origin and they often provide adaptive traits that has the ability to enhance the fitness of bacteria within a specific niche (Dobrindt et al., 2004). The putative genomic island in B. animalis subsp. lactis encodes for genes involved in cell metabolism and gene regulation and has not been found in other bacteria (Rozman et al., 2020). This could suggest that the putative genomic island including tet(W) encodes for important B. animalis subsp. lactis niche factors, which enable it to survive and compete for nutrients in the gut and has been part of the genome of B. animalis subsp. lactis long before the antibiotic era.

The tet(W), the mobile element protein and the HP gene were absent in all eight B. animalis subsp. animalis strains included in the study (Supplementary Table 1), which otherwise exhibited almost identical gene organization in the genomic region including the genes part of the putative genomic island (Figure 1). This could suggest that the tet(W), the mobile element protein and HP genes have been inserted in an ancestor of the B. animalis subsp. lactis close to subspecies differentiation and most likely lost by the three dog originating strains (2007B, 2010B, 2011B) not carrying tet(W).

The presence of a putative mobile element protein next to tet(W) has previously been reported (Ammor et al., 2008; Gueimonde et al., 2010; Rozman et al., 2020). The sequence encodes a putative DDE transposase gene that is flanked by inverted repeats upstream and downstream of 50 bp, which collectively belong to the insertion sequence (IS) 5-like element ISBian1 family that originate from B. animalis according to ISFinder (Siguier et al., 2006).

DDE transposases are able to catalyze the movement of IS elements and transposons by introducing nicks at each end of the elements (Frost et al., 2005) and are able to move within a genome or horizontally if they are part of mobile genetic element vectors such as plasmids, conjugative transposon and phages (Vandecraen et al., 2017). However, several studies have been unsuccessful in transferring tet(W) from B. animalis subsp. lactis to other species and genera (Gueimonde et al., 2010; Naghizadeh Raeisi et al., 2018; Polit et al., 2018), A BLASTp analysis showed that the IS5-like element ISBian1 family with 99.07% identity was found in the human ileum isolated Angelaksiella massiliensis (Mailhe et al., 2017) and the IS5 element was not associated with tet(W) in this species. The IS5 element was not found in other bifidobacterial species besides B. animalis subsp. lactis. The IS5 element was not found in other positions within the B. animalis subsp. lactis genomes and the inverted repeats flanking the transposase was only flanking the transposase next to tet(W). This indicates that the IS element is stably positioned next to tet(W) and does not mobilize within the B. animalis subsp. lactis genome, which is in accordance with the stable nature of the B. animalis subsp. lactis genome (Morovic et al., 2018).

Besides IS elements involvement in mobilization, IS5 elements are mainly able to modulate the expression of neighboring genes through co-transcription from the transposase promoter located in the terminal inverted repeat if inserted into non-coding regions (Schnetz and Rak, 1992; Luque et al., 2006; Vandecraen et al., 2017). The IS5 element flanking tet(W) in B. animalis subsp. lactis is positioned in a non-coding region meaning it does not cause deletion of other genes (Figure 1) and has previously been shown to be co-transcribed with tet(W) (Gueimonde et al., 2010). This indicates that the IS5 element potentially is involved in modulating the expression of tet(W) rather than mobilization.

All previously published tet(W) genes were included in the analysis. Direct submissions at NCBI also include other tet(W) genes, however, none of these exhibited 100% identity to the subspecies B. animalis subsp. lactis tet(W) and we did not find any variants not represented in the analysis (data not shown). The published tet(W) genes are therefore a good presentation of tet(W).

tet(W) is one of the most widely spread resistance genes and is both found in Gram-positive and -negative bacteria (Chopra and Roberts, 2001). Despite the wide spread nature of tet(W), it was not found to be encoded by all the strains within the examined Gram-positive and -negative species, showing that tet(W) has been acquired by a few strains or lost as compared with B. animalis subsp. lactis where it is a general genetic feature of the subspecies. For both the Gram-positive and -negative bacteria other than B. animalis subsp. lactis, tet(W) was often found to be flanked by mobile genetic elements (Table 1) and in some strains tet(W) was positioned in a genomic region with several mobile genetic elements, e.g., B. longum BG7 and A. pyogenes BBR1. Transfer of tet(W) has been reported for B. longum strain F8 (Kazimierczak et al., 2006), A. pyogenes (Billington et al., 2002), S. suis (Palmieri et al., 2011) and B. fibrosolvens (Scott et al., 1997). Within species, the tet(W) genes in the examined Gram-positive and -negative bacteria were positioned in different genomic regions. Together, this indicates that tet(W) probably has been acquired independently in the examined bacteria in Table 1.

The observation that tet(W) is generally present in B. animalis subsp. lactis strains and is positioned in the same genomic region indicates that tet(W) is conserved and thereby an innate part of the subspecies, while tet(W) only has been acquired by a few strains within the examined Gram-positive and -negative bacterial species.

A phylogenetic analysis was conducted of the tet(W) gene (Supplementary Figure 1) and protein (Figure 2) present in B. animalis subsp. lactis (Supplementary Table 1) and in the examined Gram-positive and -negative bacteria (Table 1).

The tet(W) genes encoded by the M. elsdenii strains (2–9, 7–11, 4–13) was shorter (1474–1476 bp) and exhibited a GC content (54.61–55.22%) higher compared to the other examined tet(W) genes and was therefore excluded from the phylogenetic analysis. The tet(W) gene of the remaining M. elsdenii strain (25–50) was found to be more similar to the other tet(W) genes and therefore included in the analysis.

Generally, the phylogenetic trees showed a high similarity between the different tet(W) genes and proteins, which is in agreement with previous observations (Aminov and Mackie, 2007), with the number of SNPs ranging from 1 to 46 and single amino acid polymorphisms (SAPs) ranging from 1 to 23 in the coding region compared to the tet(W) genes encoded by B. animalis subsp. lactis. The tet(W) gene encoded by A. pyogenes differed the most from B. animalis subsp. lactis tet(W) (161 SNPs and 69 SAPs). None of the SNPs lead to a premature stop codon. Based on the number of SNPs and SAPs (Table 2), clades were formed in the phylogenetic trees (Figure 2 and Supplementary Figure 1), which follows the phylogeny for B. animalis subsp. lactis but not the other examined Gram-positive and -negative bacteria.

The phylogenetic analysis showed that the tet(W) genes (Supplementary Figure 1) and proteins (Figure 2) from the B. animalis subsp. lactis strains share a high degree of homology and forms a separate clade.

The tet(W) gene and protein in the B. pseudocatenulatum type strain DSM 20438 (Genome GC content 56.40%) was located nearest the B. animalis subsp. lactis tet(W) genes and proteins in the phylogenetic trees and exhibited 12 SNPs and 5 SAPs compared to the tet(W) genes and proteins encoded by B. animalis subsp. lactis. The tet(W) gene encoded by B. pseudocatenulatum DSM 20438 and B. animalis subsp. lactis both exhibit a high identity to tet(W) from S. suis (FN396364). The tet(W) gene encoded by B. pseudocatenulatum strain 12 exhibited 45 SNPs and 22 SAPs and was located in another clade than the DSM 20438 tet(W) gene, indicating that the tet(W) encoded by the two B. pseudocatenulatum strains differ. tet(W) has been shown to be present in 33–41% of B. pseudocatenulatum isolates from human (Aires et al., 2007; Wang et al., 2017), no mobile genetic elements including IS5 elements was found in the flanking regions of tet(W) in the two examined strains (Table 1) and transfer of tet(W) from B. pseudocatenulatum have so far not been shown to occur (Wang et al., 2017). An examination of the flanking sequences of tet(W) in B. pseudocatenulatum type strain DSM 20438 revealed that the downstream genes were organized similarly as the genes downstream of tet(W) in the majority of the studied B. animalis subsp. lactis strains (Figure 1), except that six hypothetical proteins was present between tet(W) and the GMP synthase gene and no IS5-like element was present (Supplementary Figure 2). These genes were also present in B. pseudocatenulatum strain 12 but in another genetic location than tet(W), and in a B. pseudocatenulatum strain (ca_0067, NZ_RCXS00000000) that did not encode tet(W). This indicates that the presence of these genes is independent of the presence of tet(W) and are shared genes between B. animalis subsp. lactis and B. pseudocatenulatum.

The tet(W) genes present in the examined Gram-positive and -negative bacteria including the two B. pseudocatenulatum strains, were scattered over different clades in the phylogenetic tree indicating that the tet(W) genes encoded by these bacteria are diverse, does not follow the phylogeny of the specific species and thereby support the acquired nature of these tet(W) genes.

A core genome phylogenetic analysis was conducted with the examined B. animalis subsp. lactis strains (Supplementary Table 1), the bifidobacterial species from Table 1 and B. animalis subsp. animalis strains from Supplementary Table 1 (Figure 3). For each species, strains were included that both did and did not encode tet(W), except for B. animalis subsp. animalis and B. bifidum.

The core genome phylogenetic analysis showed that the bifidobacterial species separated from each other in individual clades and both strains with and without tet(W) clustered together within species, showing that the core genome analysis was able to separate at species and subspecies level.

The fact that the tet(W) gene encoded by the examined B. animalis subsp. lactis strains formed a separate clade in the gene and protein phylogenetic analysis (Supplementary Figure 1 and Figure 2) similar to the one formed in the core genome phylogenetic tree shows that the phylogeny of tet(W) follows the phylogenetic relationship of the subspecies, indicates that tet(W) originates from an ancestral host. This is further supported by the gene being positioned in the same genomic context in the examined strains. For the other examined bifidobacterial species, the tet(W) genes does not follow the phylogeny of the species, indicating that the tet(W) gene has been acquired at different timepoints, which is in line with them being flanked by different mobile genetic elements and positioned in different genomic contexts. This indicates that tet(W) present in B. animalis subsp. lactis is distinct from tet(W) found in other bifidobacterial species and other genera.