Whole genome sequencing and function prediction of 133 gut anaerobes isolated from chicken caecum in pure cultures

1Characterisation of gut microbiota is nowadays relatively simple due to recent advances in next generation sequencing. However, though DNA sequencing is useful for monitoring dynamic changes in microbiota composition, it does not enable the understanding of biological functions of individual microbiota members. Shotgun sequencing of total DNA or RNA/cDNA sequencing can partly indicate the metabolic potential of microbial communities but is limited in addressing biological functions of individual microbiota members. Even in the cases of analysis of microbial communities with low complexity when long DNA contigs can be assembled and associated with a particular bacterium, such approaches do not allow for subsequent experimental verification of observed data due to the unavailability of pure bacterial cultures. Isolation of gut anaerobes in pure cultures is therefore the best way to examine the characteristics of individual microbiota members experimentally.

2Isolation of bacterial gut microbiota members in pure cultures is usually an issue since the majority of bacterial isolates colonising the intestinal tract are strict anaerobes. Although it is commonly reported that between 10 and 50% of bacterial species colonising the intestinal tract can be grown in vitro, recent studies proposed that up to 90% of major gut colonisers may be grown in vitro if multiple culture conditions are tested. Despite this, the isolation of a particular gut anaerobe may remain an issue since even the most abundant microbiota members at species level only rarely form more than 1% of the total population. This means that the desired bacterial species may be represented by a single colony growing on an agar plate together with hundreds or thousands of others and the likelihood of picking up the particular species is indeed rather low.

3Chickens represent a specific case for studies focused on host - microbiota interactions. Chickens in commercial production are hatched in a clean environment of hatcheries without any contact with adult hens and their colonisation is dependent on environmental sources. Moreover, newly hatched chickens can be colonised by microbiota of wide a range of compositions and colonisation of the intestinal tract of commercially hatched chickens may therefore differ from the colonisation of intestinal tract of chickens which would hatch in a nest. Perhaps not surprisingly, commercially hatched chickens are highly sensitive to colonisation with different pathogens, e.g. Salmonella, nevertheless, their resistance can be increased by providing them with a complex microbiota of adult hens.

4In our recent studies, we characterised the development of gut microbiota in commercially hatched chickens throughout their whole life, identified proteins expressed by the main microbiota members and verified that gut microbiota may protect chickens against Salmonella Enteritidis infection. However, which bacterial species are the protective ones is not known and more detailed knowledge of individual microbiota members is needed. One way forward is to obtain well-characterised pure cultures of gut anaerobes. In this study we therefore cultured hundreds of isolates of anaerobes from chicken caecum and sequenced more than 200 of them. Analysis of their genomic sequences showed that we collected isolates from 7 different phyla. The aims of the subsequent comparisons focused mainly on the representatives from phyla Bacteroidetes and Firmicutes was to reveal differences in (poly)saccharide utilisation, propionate and butyrate fermentation and interactions with the host including motility or the ability to degrade host derived proteins.

5Altogether 204 isolates were obtained in pure culture and sequenced. Since in several cases we sequenced isolates exhibiting higher than 99% sequence similarity between their genomes, the final number of different isolates decreased to 133 (Table 1, Fig. 1). The lowest sequencing coverage was 43 fold for Drancourtella massiliensis An12, the highest 1526 fold for Lactobacillus gasseri An197, and median coverage over all sequenced isolates was 312 fold (see Additional file 1 for all details). The isolates belonged to 7 different phyla – Firmicutes (84 isolates), Bacteroidetes (29 isolates), Actinobacteria (15 isolates), Proteobacteria (1 isolate each of Escherichia and Desulfovibrio), Verrucomicrobia (1 isolate of Akkermansia), Elusimicrobia (1 isolate of Elusimicrobium) and Synergistetes (1 isolate of Cloacibacillus). The similarity of whole sequences of 16S rRNA to the to the most homologous GenBank entries was lower than 94% for 15 isolates. Considering taxonomic recommendations, these isolates represent candidates for new genera and two Muribaculum-like isolates might belong to a novel bacterial family or even an order. Alternative analysis based on alignment of RpoB amino acid sequences yielded similar clustering of individual isolates as that achieved by 16S rRNA comparison (Additional file 2).

6When we compared the sequence of 16S rRNA genes of all 133 isolates with the operational taxonomic units (OTUs) combined from our previous studies, rather unexpectedly, 7 isolates were excluded from the analysis by QIIME at the chimera removal step. Three isolates formed OTUs which we did not detect among the microbiota of the two studies. The rest of the isolates were assigned into particular OTUs. Nineteen isolates were assigned to 11 OTUs from the 100 most frequent OTUs, and out of the 1000 most common OTUs, we obtained 76 isolates belonging to 42 OTUs. Exact ranking of each individual isolate among the most frequent OTUs present in the chicken caecum can be found in Additional file 1.

7Genome sizes ranged from 1.51 Mb in Elusimicrobium to 6.70 Mb in Bacteroides ovatus. Larger genomes were usually recorded in isolates from phylum Bacteroidetes (Additional file 1). Actinobacteria possessed small genomes ranging from 2.1 to 2.5 Mb. Genomes of Firmicutes ranged mostly from 3 to 4 Mb although genomes of Enterococcus, [Eubacterium] cylindroides and Lactobacillus were among the smallest ones with genome sizes around 2 Mb. Genomic GC content ranged from 28.0 to 62.1% in Firmicutes, from 41.6 to 61.9% in Bacteroidetes and from 62.1 to 68.7% in Actinobacteria. GC content of individual isolates belonging to the remaining 4 phyla ranged from 50.5 to 57.9% (Additional file 1).

8Network analysis based on the correlation of individual gene counts in individual isolates confirmed similarities observed by 16S rRNA gene alignment. Individual isolates from phyla Proteobacteria, Verrucomicrobia, Elusimicrobia and Synergistetes formed disconnected vertices of the network (Fig. 2). Firmicutes split into families Lactobacillaceae, Enterococcaceae, Veillonellaceae (genera Megasphaera and Megamonas) and order Clostridiales (families Erysipelotrichaceae, Ruminococcaceae and Lachnospiraceae) except for Clostridium perfringens. Members of the family Erysipelotrichaceae formed a slightly eccentric cluster at the periphery of other isolates from the order Clostridiales indicating their slightly different coding capacity when compared to isolates belonging to families Lachnospiraceae and Ruminococcaceae. All Bacteroidetes formed a disconnected network cluster showing similar genetic composition in all isolates. The majority of Actinobacteria formed a single network cluster except for the representatives of genus Gordonibacter showing that the genome of this genus differed from the rest of Actinobacteria sequenced in this study.

9Since gut microbiota is formed mainly by representatives of phyla Bacteroidetes and Firmicutes, we specifically focused on the comparison of genomes of isolates belonging to these two phyla. Representatives belonging to phylum Bacteroidetes (Gram negative bacteria) coded for proteins required for the biosynthesis of a Gram negative cell wall while representatives of phylum Firmicutes (Gram positive bacteria) coded for proteins required for the biosynthesis of a Gram positive cell wall. However, Megamonas and Megasphaera (family Veilonellaceae, phylum Firmicutes, i.e. Gram positive bacteria) harboured genes required for the synthesis of Gram negative cell wall type (Fig. 3). Bacteroidetes and Firmicutes differed in their mode of transport across the cell wall and cytoplasmic membrane. In Bacteroidetes, genes belonging to Ton and Tol transport systems were the most frequent whilst Firmicutes encoded ABC transporters, ECF class transporters and TRAP transporters (Fig. 3). Genes enabling sporulation were specific to Gram positive Firmicutes except for Lactobacillaceae, Enterococcaceae, Veilonellaceae and two [Eubacterium] cylindroides isolates (Fig. 3). However, there were differences in the composition and distribution of genes necessary for spore formation among the isolates with sporulation potential. Members of the family Erysipelotrichaceae (Massiliomicrobiota sp. and [Clostridium] spiroforme) did not code for stage III sporulation proteins and Faecalibacterium, Anaerofilum and Gemmiger (family Ruminococcaceae) did not code for spore proteins though the rest of the genes required for sporulation were present in their genomes (Fig. 3). None of the Bacteroidetes isolates encoded proteins required for sporulation and their oxygen tolerance could be dependent on batABCDE operon (Bacteroides aerotolerance). batABE genes were present in all representatives of Bacteroidetes and batCD were present in all Bacteroidetes except for two Muribaculum isolates. When we examined the survival rate of the anaerobes after exposure to aerobic conditions experimentally, the most sensitive isolates were single isolates of Muribaculum intestinale (phylum Bacteroidetes), [Clostridium] glycyrrhizinilyticum, [Clostridium] saccharolyticum, Anaeromassilibacillus senegalensis and Flavonifractor plautii (all phylum Firmicutes) which did not survive an hour long exposure to the air. Bacteria which did not survive 8 or 24-h-long exposure to the air belonged mainly to the order Clostridiales. Within Clostridiales, 45% of the tested isolates did not survive 8-h long exposure to the air, and an additional 25% died between 8 to 24-h exposure. Only 25% of tested isolates from the order Clostridiales survived 24-h air exposure. On the other hand, representatives of Bacteroidetes and Actinobacteria were usually tolerant to sudden air exposure since 62 and 73% of tested isolates survived 24-h air exposure, respectively (Fig. 4 and Additional file 3). Bacteroides sp. encoded a high number of proteins involved in polysaccharide and monosacharide metabolism while the number of genes required for metabolism of di- and oligosaccharides was similar in both Bacteroidetes and Firmicutes (Fig. 3).

10Short chain fatty acids and butyrate in particular, are acknowledged as important metabolites of bacterial origin. All genes required for butyrate production from pyruvate and acetyl-CoA were present in the genomes of Ruminococcaceae (genera Butyricicoccus, Pseudoflavonifractor, Flavonifractor, Anaeromassilibacillus, Anaerotruncus and Faecalibacterium) (Fig. 1). In addition, this pathway was also present in Megasphaera elsdenii, [Eubacterium] cylindroides, [Clostridium] lactatifermentans, [Clostridium] saccharolyticum and [Eubacterium] hallii, the latter three belonging to the family Lachnospiraceae (Fig. 1). Butyricimonas paravirosa was the only isolate from phylum Bacteroidetes coding for all genes required for butyrate production from pyruvate and acetyl-CoA. Genes coding for enzymes of complete lysine fermentation pathway leading to butyrate production were present in three genera of the phylum Bacteroidetes (Muribaculum, Butyricimonas and Odoribacter) and the majority of Flavonifractor isolates belonging to the phylum Firmicutes. Butyricimonas and Odoribacter also encoded the whole pathway required for the conversion of succinate and 4-hydroxybutyrate into butyrate (Fig. 1). Terminal steps in butyrate production were dependent on transferases transferring CoA moiety from butyryl-CoA to acetate, acetoacetate or 4-hydroxybutyrate, or phosphate butyryltransferase (PBT) and butyrate kinase (BUK). Surprisingly, we did not find butyrate kinase in all isolates using the PBT - BUK pathway for butyrate release from butyryl-CoA. In such a case, butyrate-phosphate may serve for substrate phosphorylation in enzymatic reactions similar to acetate-phosphate.

11Propionate production via a succinate-methylmalonate pathway was characteristic of Bacteroidetes (Fig. 1) as genes for methylmalonyl-CoA mutase, epimerase and decarboxylase were detected in genomes of all isolates from this phylum. This pathway was quite rare in Firmicutes since methylmalonyl-CoA mutase, epimerase and decarboxylase were encoded only by Megamonas and two Flavonifractor isolates (Fig. 1).

12Fermentation of carbohydrates results in the production of H2 which has to be removed from the community since its increased concentration suppresses glycolysis. H2 can be removed by methanogens, acetogens or sulphate reducing bacteria. We did not isolate a single methanogen. Desulfovibrio (phylum Proteobacteria) encoded key genes for sulphate reduction to H2S (adenylylsulphate reductase, sulphate adenylyltransferase, dissimilatory sulphite reductase and sulphite reduction-associated complex DsrMKJOP). Potential for H2 removal by acetogenesis was recorded in [Clostridium] saccharolyticum, [Eubacterium] fissicatena and all Blautia isolates (B. coccoides, B. producta, B. schinkii) since all these bacteria encoded corrinoid iron-sulfur acetyl-CoA synthase and 5-methyltetrahydrofolate methyltransferase (Fig. 1).

13Since gut microbiota may interact with its host, finally we searched for the presence of genes which may facilitate such interactions (Additional file 4). Genes encoding collagenase precursor, hemagglutinin or hemolysin A were present in all 29 Bacteroidetes isolates but not in Firmicutes (Fig. 5). Hyaluronidase was present in 16 isolates from the phylum Bacteroidetes and five Firmicutes. Different heparinases were detected in 14 isolates from the phylum Bacteroidetes and seven representatives of Firmicutes – out of these all five Faecalibacterium isolates encoded heparinase II/III-like and outside this genus, heparinases were present only in two isolates from the phylum Firmicutes (i.e. Gemmiger formicilis and [Clostridium] glycyrrhizinilyticum). Chondroitinase was present in nine isolates from phylum Bacteroidetes (genera Alistipes and Bacteroides) but in no isolate from Firmicutes. Mucin-desulfating sulfatase was present in the genomes of 20 Bacteroidetes isolates (mainly in genera Alistipes, Parabacteroides and Bacteroides) but was absent from the genomes of Firmicutes. A gene for endothelin-converting enzyme 1 precursor was present in the genomes of all Bacteroidetes but not in Firmicutes. Except for two Odoribacter isolates and Parabacteroides distatonis, all the remaining representatives of the phylum Bacteroidetes encoded 3-oxo-5-α-steroid 4-dehydrogenase capable of modification of bile or steroid hormones. This gene was not detected in Firmicutes. Glutamate decarboxylase catalysing production of γ-aminobutyrate (GABA) and glutamate/GABA antiporter was present in 20 different Bacteroidetes isolates but only in two Firmicutes and these were both pathogenic Clostridium perfringens isolates. Histidine decarboxylase catalysing production of histamine was quite rare and was present only in two Bacteroides dorei isolates (Fig. 5). On the other hand, except for all Lactobacilli, the gene for fibronectin/fibrinogen-binding protein was present in all isolates from Firmicutes but in none representative from the phylum Bacteroidetes. Genes for flagellar motility were quite rare and were absent in all Bacteroidetes isolates. Complete flagellar operon was present only in two Flavonifractor isolates, one strain of Anaeromassilibacillus senegalensis and all three isolates of Anaerotruncus colihominis (Fig. 5).

14In this study we sequenced, annotated and analysed genomes of 133 different anaerobes cultured from the chicken caecum. Several isolates represented species which are available in only a few pure cultures worldwide – a single manuscript reported a pure culture of Elusimicrobium and only two papers reported culture of Cloacibacillus porcorum. Although such isolates are of clear potential for future experiments, in this report we mainly focused on the comparison of genetic potential of isolates belonging to two main phyla inhabiting the intestinal tract of chickens, Bacteroidetes and Firmicutes.

15Whole genome sequencing showed that representatives of Bacteroidetes and Firmicutes have a genetic potential to follow different strategies of gut colonisation which may explain their coexistence in the intestinal tract. Bacteroidetes encoded genes for the biosynthesis of a Gram negative cell wall, Gram positive anaerobes encoded genes for the biosynthesis of a Gram positive cell wall. Only Megamonas and Megasphaera, though belonging to Gram positive Firmicutes, encoded genes for the biosynthesis of Gram negative cell wall type, in agreement with previous report. The Ton/ExbD transport system of Bacteroidetes has been identified earlier as highly expressed in vivo and ABC, ECF and TRAP transporters were described as characteristic of Firmicutes. Bacteroidetes were reported to increase in gut microbiota with high fiber content in their diet and to forage on host derived polysaccharides in the absence of fiber. Consistent with this, Bacteroidetes and the genus Bacteroides in particular encoded a high number of genes required for polysaccharide metabolism. Since the polysaccharides of feed and host origin consist of various (amino)monosaccharides, Bacteroidetes encoded also a wider range of genes required for monosaccharide degradation than Firmicutes (Fig. 3).

16Butyrate is produced mainly by Firmicutes. Carbohydrate fermentation to pyruvate and acetyl-CoA was the most frequent butyrate production pathway, as proposed earlier. Butyrate production was mainly associated with Ruminococcaceae and less frequently with Lachnospiraceae or Erysipelotrichaceae. Bacteroidetes are not the main butyrate producers as acetyl-CoA conversion to butyrate was only found in Butyricimonas. However, Bacteroidetes were capable of butyrate production by alternative pathways, e.g. from 4-hydroxybutyrate as recorded in Butyricimonas and Odoribacter or by lysine fermentation as recorded in Muribaculum, Butyricimonas and Odoribacter. This is in line with conclusions derived from human microbiota studies.

17Gut microbiota interacts with the host. The potential of Firmicutes to interact with the chicken host seems to be less extensive in comparison to Bacteroidetes. Except for Lactobacilli, Firmicutes isolates encoded fibronectin/fibrinogen-binding protein. Interestingly, we detected chicken fibrinogen-domain containing proteins as tightly associated with gut microbiota. Such proteins aggregate bacteria and enable association of different gut microbiota members, based on current results, preferentially those belonging to the phylum Firmicutes. Binding to chicken fibrinogen-domain containing proteins may result in the formation of random bacterial aggregates and those with the most optimal composition, e.g. butyrate producers releasing H2 with Blautia consuming H2 for acetate production thus allowing glycolysis in butyrate producers to proceed, will rapidly multiply and define final microbiota composition. Except for fibrinogen binding, representatives from the phylum Bacteroidetes had a greater potential to affect host behaviour than representatives from the phylum Firmicutes. Bacteroidetes encoded collagenase, hemagglutinin, hemolysin, hyaluronidase, heparinases, chondroitinase or mucin-desulfating sulfatase required for degradation of host structures. In addition, representatives from the phylum Bacteroidetes encoded endothelin-converting enzyme 1 precursor, which plays a significant role in cardiovascular diseases and Alzheimer’s disease in humans, 3-oxo-5-α-steroid 4-dehydrogenase capable of modification of steroid hormones or bile, or glutamate decarboxylase catalysing production of γ-aminobutyrate (GABA), a mediator within the enteric nervous system. However, it will have to be elucidated whether these genes are expressed in vivo and whether their products may reach host structures.

18Finally, the rather counterintuitive conclusion came from the prediction of survival following air exposure. Although Bacteroidetes should be more sensitive to air exposure than Firmicutes which are capable of spore formation, our data clearly showed that Clostridiales, including all important butyrate producers, were highly sensitive to air exposure. Due to the experimental design, we likely determined sensitivity of vegetative cells and not of spores. Despite this, the extreme sensitivity of vegetative cells of Clostridiales may explain their commonly reported disappearance in inflammatory bowel disease patients. Inflammation leads to locally increased oxygen levels due to the oxidative burst of granulocytes and macrophages. Disappearance of Clostridiales including major butyrate producers can therefore be a consequence rather than a cause of inflammatory bowel diseases. Similarly, a reported increase in the abundance of Bacteroidetes or Megasphaera during inflammatory diseases may be a mere consequence of their higher resistance to oxygen and the disappearance of oxygen sensitive bacterial species from order Clostridiales. An overgrowth of facultative anaerobes like those from family Enterobacteriaceae in acute colitis also fits in the proposed scenario.

19In this study we isolated and sequenced 133 different strains originating from chickens intestinal tracts belonging to seven different phyla. Analysis of their genomic sequences showed that butyrate production was mainly associated with Ruminococcaceae, and less frequently with Lachnospiraceae or Erysipelotrichaceae, all belonging to phylum Firmicutes. Representatives of phylum Bacteroidetes commonly encoded proteins (collagenase, hemagglutinin, hemolysin, hyaluronidase, heparinases, chondroitinase, mucin-desulfating sulfatase or glutamate decarboxylase) that may enable them to interact with their host. Even such a brief list of genes shows that representatives of Bacteroidetes and Firmicutes follow different strategies of gut colonisation which contributes to their coexistence in the intestinal tract.

20The handling of animals in this study was performed in accordance with current Czech legislation (Animal Protection and Welfare Act No. 246/1992 Coll. of the Government of the Czech Republic). Experiments performed in this study were approved by the Ethics Committee of the Veterinary Research Institute (permit number 4/2016) followed by the Committee for Animal Welfare of the Ministry of Agriculture of the Czech Republic.

21The chickens were sacrificed under chloroform anesthesia by cervical dislocation. Whole caeca with their contents originating from 18 random healthy chickens or hens 4 to 40 weeks of age were removed during necropsy, chilled on ice and transported to an anaerobic chamber for further processing within one hour. The caeca were opened in an anaerobic chamber (10% CO2, 5% H2 and 85% N2 atmosphere; Concept 400, Baker Ruskinn, USA) and 0.5 g of content was squeezed into 4.5 ml pre-reduced PRAS dilution blank (0.1 g magnesium sulfate heptahydrate, 0.2 g monobasic potassium phosphate, 0.2 g potassium chloride, 1.15 g dibasic sodium phosphate, 3.0 g sodium chloride, 1.0 g sodium thioglycolate, 0.5 g L-cysteine, 1000 ml distilled water; final pH 7.5 +/− 0.2 at 25 °C) and mixed thoroughly. All samples were serially diluted in pre-reduced PRAS dilution blank and plated on Wilkins-Chalgren anaerobe agar (WCHA) (Oxoid) supplemented with 30% of rumen fluid. The rumen fluid was collected from cows by an oral probe, filtered through cheesecloth, centrifuged at 8000 g for 30 min and sterilised by filtration through a 0.22 μm filter. Aliquots of rumen fluid were stored at − 20 °C. WCHA was additionally supplemented with 5 mg/l hemin, 1 mg/l cellobiose, 0.5 g/l soluble starch, 1 mg/ml maltose, 0.2 ml vitamin K1 solution (0.1 ml of filter sterilized vitamin K1 in 20 ml 95% ethanol) and 0.5 mg/ml L-cysteine. Approx. 10 well-separated colonies of different morphology were selected from each agar plate after a five-day incubation at 37 °C and purified by subculture on WCHA. All isolates were stored at − 80 °C in pre-reduced PRAS dilution blank containing glycerol at 20% concentration and equal volume of sterile sheep blood. Sensitivity of pure anaerobe cultures to air oxygen exposure was tested exactly as described elsewhere. Briefly, bacterial cultures were serially diluted in anaerobic chamber and plated on 4 copies of WCHA. One copy of WCHA was left in the anaerobic chamber to determine initial counts of each anaerobe. The remaining 3 copies of WCHA plates were placed into a standard aerobic 37 °C incubator and after 1, 8 and 24 h, a single copy of agar plate was returned back to the anaerobic chamber to check for growth restoration.

22DNA was purified using DNeasy Blood & Tissue Kit (Qiagen). Sequencing library was prepared from 1 ng of RNA-free genomic DNA using the Nextera XT DNA Sample Preparation kit (Illumina) and whole genome sequencing was performed using the NextSeq 500/550 High Output Kit v2 and Illumina NextSeq 500 sequencing platform in the paired-end modus (2 × 150 bp). Raw sequencing reads were quality trimmed using Trimmomatic v0.32 with the sliding window of 4 bp and average quality threshold value equal to 17. Minimal read length was set to 48 bp. Trimmed paired-end reads were assembled via de novo assembler IDBA-UD v1.1.1 with k-mer sizes ranging from 20 to 110 with an increment of 15. Contigs with coverage lower than 10% of average coverage of L50 contigs were filtered out and the remaining contigs were scaffolded employing SSPACE scaffolder v3.0. All scaffolds containing N’s in their sequences were split into N-free sequences. Finally, scaffolds with a length shorter than 500 bp were discarded.

23Species definitions used in this study are based on the BLAST comparison of whole 16S rRNA sequences against entries deposited in NCBI 16S rRNA sequence database performed on January 4, 2017. For clarity of the paper, we used the designation of the most similar bacterium based on the lowest E-value for description of our isolates, thus apparently ignoring the fact that in some cases there was 100% identity whilst in the opposite extreme, the sequence of 16S rRNA of one of our isolates was only 83% similar to the closest relative deposited in the NCBI 16S rRNA database (Additional file 1). All 16S rRNA sequences were compared also against RDP SeqMatch database (on January 10, 2017) which allowed for alternative taxonomy including classification of individual isolates into higher taxonomic units. In addition, ribosomal protein multilocus sequence typing (rMLST) and GTDB organism identification (http://gtdb.ecogenomic.org/) databases were used to verify taxonomic classification. Taxonomic classification by these alternative protocols is provided in Additional file 1. Gene predictions and functional annotations were performed by RAST. Assembled and annotated genomes as well as raw sequencing data were deposited in NCBI under accession number PRJNA377666 and genomes with comprehensive RAST annotation are available upon request.

24To exclude genomes of the same isolates picked up on independent occasions from the subsequent analyses, whole genome sequences of all isolates were mutually compared using the QUAST tool v3.1. Two genomes were considered as identical if they shared more than 99% of genome content and had less than 1 indel per 100 kb. A single isolate was selected as a representative of each group of identical isolates for all downstream analyses.

25Whole gene content similarity clustering was compared using Gene Co-Expression Network analysis. This protocol detects genes with similar transcriptional regulatory program (potential members of the same pathway, protein complex, etc.). In our case, the gene expression vector for particular gene was replaced by gene copy-number vector for particular bacteria and the protocol therefore detected bacterial isolates with similar gene content. Interconnected bacteria shared more than 50% of genes based on gene name designation provided by RAST annotation. Whole gene content similarity of individual isolates was analysed in R. At first, matrix of Spearman’s correlation coefficients was calculated for all pairs of isolates using vectors of respective gene counts. The correlation matrix was then transformed to the adjacency matrix using threshold of + 0.5 as a cut-off for two vertices to be considered as connected. The undirected network was constructed from such a matrix using igraph package (http://igraph.org/r/), and edge betweenness based community structure detection algorithm was then employed to identify separate network modules. Communities with more than three members were considered nontrivial and were highlighted in a network plot.

26Trimmed reads were aligned against SILVA bacterial 16S rRNA database using SortMeRNA v2.1 and extracted 16S rRNA reads were assembled via de Bruijn graph-based de novo assembler SPAdes v3.6.0. Finally, sequences coding for 16S rRNA genes were predicted employing barrnap tool v0.6 (https://github.com/tseemann/barrnap). For the purposes of phylogenetic analysis, 16S rRNA genes were aligned by ClustalW v2.1 using default gap penalties, DNA weight matrix IUB and transition weight 0.2. The phylogenetic tree topologies were inferred employing Bayesian statistics via MrBayes v3.2.6 using the parameters as follows: mixed model of nucleotide substitution, gamma-distributed rates among sites, four Monte Carlo Markov chains for 7,000,000 cycles which were sampled every 1000th generation and the first 25% of the samples were discarded as burn-in. The final tree topology was generated employing 50% majority-rule consensus. Average standard deviation of split frequencies was 0.0083, maximum standard deviation of split frequencies was 0.1095, average potential scale reduction factor was 1.000 and maximum potential scale reduction factor was 1.015. The final tree topology was visualized by iTOL v3.4.3 and edited using Inkscape v0.91 (www.inkscape.org).