Bacterial hosts of clinically significant beta-lactamase genes in Croatian wastewaters

Abstract

Wastewater treatment plants (WWTPs) provide a suitable environment for the interaction of antibiotic resistant bacteria and antibiotic-resistance genes (ARGs) from human, animal, and environmental sources. The aim was to study the influent and effluent of two WWTPs in Croatia to identify bacterial hosts of clinically important beta-lactamase genes (bla_TEM, bla_VIM, bla_OXA-48-like) and observe how their composition changes during the treatment process. A culture-independent epicPCR (Emulsion, Paired isolation and Concatenation Polymerase Chain Reaction) was used to identify the ARG hosts, and 16S rRNA amplicon sequencing to study the entire bacterial community. Different wastewater sources contributed to the significant differences in bacterial composition of the wastewater between the two WWTPs studied. A total of 167 genera were detected by epicPCR, with the Arcobacter genus, in which all ARGs studied were present, dominating in both WWTPs. In addition, the clinically important genera Acinetobacter and Aeromonas contained all ARGs examined. The bla_OXA-48-like gene had the highest number of hosts, followed by bla_VIM, while bla_TEM had the narrowest host range. Based on 16S rRNA gene sequencing, ARG hosts were detected in both abundant and rare taxa. The number of hosts carrying investigated ARGs was reduced by wastewater treatment. EpicPCR provided valuable insights into the bacterial hosts of horizontally transmissible beta-lactamase genes in Croatian wastewater.

16S rRNA, antibiotic-resistance genes, bacteria, Croatia, epicPCR, wastewater

Introduction

Antimicrobial resistance (AMR) is one of the major health issues of the modern era. World Health Organization has ranked AMR as one of the top 10 threats for the global health (WHO 2019). Wastewater treatment plants (WWTPs) are of particular interest for the study of AMR because they collect wastewater from multiple sources, including human and non-human sources, and provide a suitable environment for close bacterial interaction and exchange of antibiotic-resistance genes (ARGs). In WWTPs, high nutrient concentration, high bacterial density, and the presence of sub-inhibitory concentrations of antibiotics and other selective agents such as heavy metals and disinfectants favour horizontal gene transfer (HGT) in bacteria (Rizzo et al. 2013, Karkman et al. 2018, Manaia et al. 2018, Pazda et al. 2019). Of particular interest is to know which bacteria are the hosts of clinically significant ARGs and how the composition of the hosts is altered by the wastewater treatment process. EpicPCR (Emulsion, Paired Isolation and Concatenation Polymerase Chain Reaction) allows us to link a gene of interest to its host bacterium in a culture-independent way (Spencer et al. 2016). Previously, this technique has been used to determine the bacterial host diversity of the dissimilatory sulfite reductase gene dsrB from lake water (Spencer et al. 2016), two ARGs (tetM, bla_OXA-58), and the class 1 integron-integrase gene (intI1) from wastewater (Hultman et al. 2018), or the integrative and conjugative XT/R391 elements from river water (Roman et al. 2021). All of these previous studies have demonstrated the efficiency of epicPCR in the analysis of complex environmental samples.

The aim of this study was to investigate by epicPCR the hosts of the beta-lactamase genes bla_TEM, bla_VIM, and bla_OXA-48-like from the WWTPs of a large (Zagreb) and a small (Varazdin) continental city in Croatia. These genes were selected based on their high clinical significance and poor removal during wastewater treatment, as determined in a previous study by Puljko et al. (2022) using quantitative PCR (qPCR). Beta-lactamase bla_TEM genes contain variants belonging to the extended-spectrum beta-lactamases that confer resistance to antibiotics of importance in human medicine, such as third generation cephalosporins (Castanheira et al. 2021). bla_OXA-48-like and bla_VIM are carbapenemase genes that hydrolyse almost all beta-lactams, including carbapenems, which are used as a last resort in the treatment of infections caused by multidrug-resistant Gram-negative pathogens (Queenan and Bush 2007). All of the above genes are commonly found on plasmids, making them easily transferable between environmental bacteria and human and animal pathogens that could be induced in WWTPs.

Materials and methods

Sampling

Two WWTPs in two continental cities (Zagreb and Varazdin) in Croatia were investigated. Both WWTPs have primary treatment, where larger particles are removed mechanically, and secondary treatment with activated sludge. The Zagreb WWTP is the largest WWTP in the country. It is located in the capital and is designed for 1 200 000 population equivalents and receives domestic and hospital wastewater (nine hospitals). The Varazdin WWTP is designed for 140 000 population equivalents and receives domestic, hospital, and industrial wastewater (dairy and poultry industry). The 24-h composite samples of the influent and effluent were taken by the technical staff of the two WWTPs on three consecutive days (28 September 2021, 29 September 2021, and 30 September 2021). These samples were considered as biological replicates. Samples were collected in sterile 500 ml (influent) and 2.5 l (effluent) bottles and brought to the laboratory within 3 h of collection. Samples were then immediately filtered and centrifuged as described below and further analysed using methodologies outlined in Fig. 1.

Figure 1.

Outline of the research methods used in this study.

DNA extraction and 16S rRNA amplicon sequencing

For DNA extractions, 50–60 ml of influent and 100–400 ml of effluent were filtered in triplicate through sterile membrane filters with a pore size of 0.2 µm (ME 24/21 ST, diameter 47 mm, Whatman, GE Healthcare, Life Science, USA). The filters were cut into smaller pieces with sterile scissors and placed in tubes for DNA extraction using the PowerSoil Kit (Qiagen, USA). DNA quantity was determined using Qubit Fluorometer 3.0 (Thermo Fisher Scientific, USA), and quality (ratio 260/280) was determined using Nanodrop Spectrophotometer (BioSpec Nano, Shimadzu, Japan). The DNA extracts were stored at −20°C.

PCR amplification of the V3–V4 region was performed with primers 341F1-4 and 785R1-4 containing Illumina TruSeq adapter sequences at the 5′ ends using Phusion Hot Start Polymerase with GC buffer (New England Biolabs, MA, USA). PCR conditions and a list of primers are given in the Supplementary material (Table S1). Amplicons were sequenced on the Illumina MiSeq platform at the Institute of Biotechnology, University of Helsinki, Finland. All amplicon sequences have been submitted to the Sequence Read Archive (PRJNA1012856).

EpicPCR

To collect cells for epicPCR, 2 ml of the influent was centrifuged at 11 000 × g for 5 min, and 200 ml of the effluent was centrifuged at 8000 × g for 40 min. Cell pellets were stored in 20% glycerol at −20°C. The ARGs selected in this study are bla_TEM, bla_VIM, and bla_OXA-48-like. Prior to performing epicPCR, conventional PCR was performed using the primers and conditions described previously (Puljko et al. 2023) to confirm the presence of the tested genes in all extracted DNA from the collected samples. At the Zagreb WWTP, the presence of bla_VIM and bla_TEM was confirmed in the influents and effluents. However, the presence of bla_OXA-48-like was confirmed by PCR in the influent, but not in the effluent. At the Varazdin WWTP, bla_TEM and bla_OXA-48-like were confirmed in the influent and effluent, but bla_VIM only in the effluent. Only the genes that were detected in both the influent and effluent were analysed with epicPCR.

EpicPCR was performed according to the original protocol described by Spencer et al. (2016) with modifications described by Hultman et al. (2018) and Roman et al. (2021). To obtain an estimate of the number of cells to be used for the epicPCR protocol (1−2 × 10⁷ cells), cells were stained with SybrGreenII (Thermo Fisher Scientific, Waltham, USA) and counted using a Petroff-Hausser counting chamber and a fluorescence microscope (Nikon eclipse e600, Japan). After counting, cells were encapsulated into polyacrylamide beads as described in Hultman et al. (2018). The presence of single cells in individual beads was confirmed by fluorescence microscopy. Fusion PCR was then performed to amplify and merge the desired ARG and a part of the 16S rRNA gene as described in Spencer et al. (2016). For fusion PCR, Phusion Hot Start Polymerase with GC buffer (New England Biolabs, USA) was used. After purification of PCR fusion products using a Monarch PCR and DNA Cleanup Kit (New England Biolabs, USA), blocking PCR was performed as described by Roman et al. (2021) to remove non-fused products. The blocking PCR products were purified and a final nested PCR was performed using Illumina TruSeq adapters to increase specificity within the fused products (Spencer et al. 2016). The nested PCR products were sequenced on the Illumina MiSeq platform at the Institute of Biotechnology, University of Helsinki, Finland. The list of all primers and PCR conditions is provided in the Supplementary material (Table S2). All epicPCR sequences have been submitted to the Sequence Read Archive (PRJNA1012856).

16S rRNA gene data analysis

Cutadapt v.3.5 (Martin 2011) was used to remove primers from the 16S rRNA sequences. FastQC v.0.11.8 (Andrews 2010) and multiqc v.1.12 (Ewels et al. 2016) were used for quality checking. Amplicon sequence variant (ASV) inference was done in R Studio v.4.2.2 (Posit team 2023) with DADA2 v.1.26.0 (Callahan et al. 2016) pipeline. Trimming length was set to 280 bp for forward reads and 220 bp for reverse reads, and maximum expected errors were set to 2 for both reads. Chimeras were removed using removeBimeraDenovo and taxonomic classification against the Silva database v.138.1 (Quast et al. 2013) was performed using assignTaxonomy function. ASVs annotated as Mitochondria and Chloroplasts were removed from the data before any further analysis.

All statistical analyses were performed in R Studio v.4.2.2 (Posit team 2023). Alpha diversity was assessed by calculating the Shannon diversity index using the phyloseq package v.1.42.0 (McMurdie and Holmes 2013). The comparison between the alpha diversity of the WWTPs was carried out using the Wilcoxon rank sum test. To compare the alpha diversity of different sample types using Shannon diversity index values, a linear mixed effects model was created using the lme4 v1.1.33 package (Bates et al. 2015). Sample type was considered a fixed effect, and WWTPs were considered a random effect. An analysis of variance was performed to assess if the constructed model was statistically significant. To assess beta diversity, samples were transformed to proportions. Principal coordinate analysis (PCoA) plot using the Bray–Curtis dissimilarity index was generated using phyloseq v.1.42.0 (McMurdie and Holmes 2013). Comparisons of bacterial composition between samples were performed with permutational multivariate analysis of variance (PERMANOVA) and the Bray–Curtis dissimilarity index using the adonis2 function from the vegan package v.2.6.4 (Oksanen et al. 2022).

EpicPCR data analysis

The quality of epicPCR sequences was checked using FastQC v.0.11.8 (Andrews 2010) and multiqc v.1.12 (Ewels et al. 2016). Demultiplexing, primer removal, quality filtering, and sequence spliting was performed using cutadapt v.3.5 (Martin 2011). Minimum length of the reads was set to 150 bp, overlap to 10 bp and error rate to 0.25. The quality cutoff values for forward reads were 5 and 20 for reverse reads. Only reads containing both the 16S rRNA and ARG parts were filtered (to exclude false positives) and split into two separate fastq files based on functional gene sequences. The ARG sequences were verified with BLASTN against the NCBI database to ensure that the correct ARG was amplified. Annotation of the 16S rRNA gene sequence was performed using mothur v.1.48.0 with the classify.seqs command (Schloss et al. 2009) against the Silva database v.138.1 (Quast et al. 2013). Given the qualitative nature of epicPCR resulting from many cycles in three PCR steps, data were analysed based on the presence/absence of a gene in a taxon. A taxon was considered present if it was found in all three biological replicates. All visualizations were created using the ggplot2 v.3.4.2 package (Wickham 2016).

Results

Bacterial composition in wastewater

The composition of the bacterial community in the Zagreb WWTP was characterized by a high relative abundance of the Arcobacteraceae family (40 ± 11%) in the influent, while the Flavobacteriaceae (27 ± 7%) and Moraxellaceae (16.3 ± 0.2) families dominated the effluent (Fig. 1). The two most abundant families in the influent of the Varazdin WWTP were Arcobacteraceae (19 ± 10%) and Streptococcaceae (23 ± 8.2%), while in the effluent, the families Rhodocyclaceae (27.2 ± 1.2%) and Arcobacteracae (25 ± 3.3%) were dominant (Fig. 2). Based on the Shannon diversity index, the wastewater from the Varazdin WWTP had slightly higher alpha diversity compared to the Zagreb WWTP, but without statistical significance (Supplementary material, Fig. S1). The bacterial community of the influents of both WWTPs was more diverse than that of the effluents (linear mixed effects model, P = 0.04). The bacterial composition of the influent and effluent of the same WWTP appeared to be significantly different (P = 0.008, R² = 32.33), as well as between the two WWTPs studied (P = 0.01, R² = 31.56), based on PERMANOVA using the Bray–Curtis dissimilarity index (Fig. 3).

Figure 2.

Bacterial composition of wastewater from Zagreb and Varazdin WWTPs in influent (IN1-3) and effluent (EF1-3) based on analysis of 16S rRNA amplicon sequences. Top 20 families with the highest relative abundance are presented.

Figure 3.

Principal coordinate analysis plot of the relative abundance of bacteria in wastewater samples using the Bray–Curtis dissimilarity index. The significant difference of all sample types from each other was confirmed with permutational multivariate analysis of variance.

Hosts of the studied ARGs

A total of 167 bacterial genera were identified as potential hosts of ARGs tested by epicPCR: bla_TEM, bla_VIM, and bla_OXA-48-like. A list of all taxa detected by epicPCR is given in Supplementary material (Table S3). In Zagreb WWTP, bla_VIM was detected in more hosts in influents (103 genera) than in effluents (40 genera), while bla_TEM was detected in almost equal numbers in both wastewater types (26 influent, 28 effluent) (Fig. 4). In addition, the ARGs bla_TEM and bla_VIM shared 25 host genera in the influent and 19 genera in the effluent of the Zagreb WWTP.

Figure 4.

Number of genera detected by epicPCR carrying antibiotic-resistance genes studied. The genes bla_OXA-48-like (in the effluent of the Varazdin WWTP) and bla_VIM (in the influent of the Zagreb WWTP) had the highest number of hosts.

In the Varazdin WWTP, however, higher number of hosts carrying bla_TEM was found in the influent (61 genera) than in the effluent (21 genera) (Fig. 4). The read quality of the sequences of bla_OXA-48-like in the influent was of poor quality and was discarded for further analysis. Therefore, only the bacterial hosts of the bla_OXA-48-like gene from the Varazdin WWTP effluent were presented here, which were represented in 106 genera, making it the gene with the highest number of hosts in this study (Fig. 4). The ARGs bla_TEM and bla_OXA-48-like shared 20 host genera in the effluent of the Varazdin WWTP.

In terms of bacterial abundance, the most dominant genus carrying all the genes studied in both WWTPs was Arcobacter (Fig. 5). Furthermore, the genera Acinetobacter and Aeromonas were also abundant carriers of all investigated genes in both WWTPs. Other significant ARG-carrying genera in the Zagreb WWTP in terms of relative abundance were Alkanindiges, Bacteroides, Flavobacterium, Streptococcus, and Trichococcus (Fig. 5). In the Varazdin WWTP, additional abundant ARG hosts were Dechloromonas spp., Lactococcus spp., Pseudarcobacter spp., Streptococcus spp., and Trichococcus spp. (Fig. 5).

Figure 5.

Bacterial genera detected by 16S rRNA amplicon sequencing. The size of the circle indicates relative abundance. Only genera with a relative abundance of ≥0.5% are shown. If the investigated antibiotic-resistance gene was detected in a particular genus by epicPCR, it is marked with a full circle, otherwise the circle is empty.

The genera detected by epicPCR that contain clinically significant pathogens are shown in Table 1. Interestingly, all clinically significant genera listed in Table 1 were the hosts of bla_OXA-48-like. In addition, there were some genera that were present in both influent and effluent (based on 16S rRNA gene sequencing) but carried studied ARG only in the effluent (based on epicPCR), which could indicate possible HGT events. In the Zagreb WWTP, such events were detected for both bla_VIM and bla_TEM in 7 and 10 genera, respectively. In the Varazdin WWTP, four events were recorded for bla_TEM. There are no data on possible HGT events of bla_OXA-48-like, as the hosts were not detected in the influent with epicPCR.

Table 1.

Clinically significant bacterial genera detected with epicPCR that carry studied antibiotic-resistance genes.

Clinically significant genera detected with epicPCR	Zagreb WWTP				Varazdin WWTP
	blaVIM		blaTEM		blaTEM		blaOXA-48-like
	IN	EF	IN	EF	IN	EF	EF
Acinetobacter	+	+	+	+	+	+	+
Aeromonas	+	+	+	+	+	+	+
Arcobacter	+	+	+	+	+	+	+
Bacteroides	+		+		+	+	+
Clostridium sensu stricto 1		+			+		+
Enterobacter							+
Escherichia-Shigella	+	+	+	+	+	+	+
Flavobacterium	+	+			+		+
Klebsiella	+	+	+		+		+
Prevotella							+
Pseudomonas	+	+	+	+	+		+
Shewanella	+					+	+
Streptococcus	+	+	+	+	+		+
Yersinia							+

Discussion

In this study, the bacterial composition of wastewater from two Croatian WWTPs (influent and effluent) was investigated using 16S rRNA amplicon sequencing. In parallel, epicPCR was used to identify host bacteria of clinically important beta-lactamase genes. The bacterial composition of the influent wastewater changed during the treatment process in both WWTPs, as evidenced by significant differences between the influent and effluent wastewater in terms of alpha and beta diversity. The bacterial community of the raw influent wastewater was more diverse (higher Shannon diversity index values) than that of the treated effluent, which is contrary to several previous reports (Hultman et al. 2018, Azli et al. 2022). The reasons for this opposite observation could be different physical and chemical variables that occur during the wastewater treatment process, resulting in lower species diversity and abundance (Do et al. 2019). In addition, the PCoA plot and PERMANOVA tests revealed that the influent was significantly different from that of the effluent at each WWTP and between WWTPs. The different composition of wastewater samples (both influent and effluent) between the two WWTPs was to be expected given the different size of the cities and the influence of industrial wastewater in the Varazdin WWTP.

The 16S rRNA gene analysis further revealed that the Arcobacteraceae family, especially Arcobacter spp. dominated in the Zagreb WWTP influent. Hultman et al. (2018) recorded the same dominance of Arcobacteraceae in raw influents of Finnish WWTPs. Members of the family Arcobacteraceae are found in diverse aquatic environments with the highest prevalence in sewage (Venâncio et al. 2022). Among the genus Arcobacter, there are several clinically important species that cause infections in humans, such as A. butzleri and A. cryaerophilus (Venâncio et al. 2022), making the detection of ARGs in Arcobacter particularly interesting. Moreover, in this study, all three targeted beta-lactamase genes (bla_VIM, bla_OXA-48, and bla_TEM) were detected by epicPCR in Arcobacter spp., suggesting that this genus may serve as an unexplored reservoir for clinically important ARGs. Based on the 16S rRNA gene sequencing results presented here, the wastewater treatment process at the Zagreb WWTP appears to effectively reduce the relative abundance of bacteria from the Arcobacteraceae family. In contrast, at the Varazdin WWTP, the abundance of Arcobacteraceae in the treated wastewater slightly increased, indicating the different environmental conditions favouring the survival of this family by the wastewater treatment process.

Flavobacteriaceae and Moraxellaceae were also very abundant families in the effluent of the Zagreb WWTP. Both families contain members that are typical environmental bacteria found in various environments, as well as opportunistic pathogens for humans and animals. Members of the genus Flavobacterium contained bla_VIM and bla_OXA-48-like genes in the effluents of the two WWTPs studied. In the Moraxellaceae family, the genus Acinetobacter, which carried all beta-lactamase genes studied, was particularly abundant. Some members of Acinetobacter spp. such as A. junii are important populations in the wastewater treatment process because they can effectively accumulate polyphosphates (Han et al. 2018). On the other hand, A. baumannii is a clinically important emerging opportunistic pathogen in humans with a high level of antibiotic resistance, and multidrug resistant strains were previously isolated from the Zagreb WWTP by Hrenovic et al. (2016). Another significant genus that contained all of the ARGs examined (bla_VIM, bla_OXA-48, and bla_TEM), as shown by the epicPCR in this study, is Aeromonas. Members of this genus are ubiquitous in aquatic environments and thrive in water distribution systems due to their ability to form strong biofilms (Batra et al. 2016). In addition, the species A. hydrophila, A. caviae, A. veronnii, and A. sobria are emerging human pathogens that cause gastrointestinal infections (Batra et al. 2016, Bello-López et al. 2019). Aeromonas spp. are already known to possess a variety of ARGs localized on mobile genetic elements and they readily exchange genes with genetically distant bacteria (Bello-López et al. 2019). In a previous study by Drk et al. (2023), the plasmid-localized carbapenemase genes bla_VIM and bla_OXA-48-like were detected for the first time in Aeromonas isolates, particularly pathogenic A. caviae, isolated from effluent samples of the Zagreb WWTP. Together with the epicPCR detection of these genes in Aeromonas in this study, this confirms that Aeromonas spp. with carbapenem-resistant genotypes remaining in wastewater after treatment could pose a potential public health threat if transmitted to humans via wastewater.

The family Streptococcacae, in particular the genera Lactococcus and Streptococcus (group of lactic acid bacteria) were among the most abundant species in the influent of the Varazdin WWTP and probably originated from the wastewater of the neighbouring dairy industry. In this work, members of both genera were found to carry bla_TEM and bla_OXA-48-like genes by epicPCR, while bla_VIM was only detected in Streptococcus spp. In addition to Arcobacteraceae, effluent from the Varazdin WWTP was dominated by Rhodocyclaceae, a family of environmental bacteria with diverse lifestyles. The most abundant genus, Dechloromonas, is common in WWTPs worldwide and is considered an important genus capable of degrading aromatic compounds (Petriglieri et al. 2021). However, this work showed that Dechloromonas spp. carried all of the ARGs studied (bla_VIM, bla_OXA-48-like, and bla_TEM), suggesting that this genus may serve as a reservoir for the dissemination of horizontally transmissible ARGs of clinical importance in the environment.

The number of bla_VIM and bla_TEM gene hosts decreased by the wastewater treatment process, which is consistent with Hultman et al. (2018). Based on the quantification of ARGs by qPCR in wastewater from the same WWTPs, Puljko et al. (2022) in 2020 detected no statistically significant difference in the relative abundance of targeted genes between influent and effluent. Therefore, the decline in the host range of bla_VIM and bla_TEM is probably due to the reduced abundance of hosts carrying the ARGs. Based on the 16S rRNA gene detection, Puljko et al. (2022) reported a 1.54 ± 0.06 (in Zagreb) and 0.71 ± 0.1 (in Varazdin) log reduction of total bacteria in wastewater in the summer after passing through the wastewater treatment process. Additionally, Puljko et al. (2022) observed lower relative abundance of bla_VIM (in Zagreb) and bla_OXA-48-like (in Varazdin) by qPCR compared to bla_TEM. However, in the present study, bla_VIM and bla_OXA-48-like were detected in more genera than bla_TEM, indicating their broader host range. This suggests that higher gene relative abundance does not imply broader host range. Possible HGT events were identified for bla_VIM and bla_TEM. However, it is unclear if all those events are indeed HGT or if some genera are very rare and therefore below detection limit of the method.

The bla_OXA-48-like gene is a highly mobile carbapenemase gene commonly found in the Enterobacteraceae family (Boyd et al. 2022). However, in this study, bla_OXA-48-like was detected in 54 different families, including potential pathogens (e.g. Acinetobacter spp., Klebsiella spp., Escherichia–Shigella spp., Enterobacter spp.), indicating a broad host range in wastewater and highlighting the potential of WWTP effluent as an environment for the spread of bla_OXA-48-like genes through vertical and horizontal transmission. This gene was also recently found in enterobacterial isolates belonging to clinically important genera of the family Enterobacteraceae, such as Escherichia and Klebsiella, in effluent water from the Zagreb WWTP (Puljko et al. 2024).

Due to the current limitation of epicPCR, it is not possible to determine the ARG hosts down to the species level because the small V4 region of the 16S rRNA gene is amplified. In addition, the method only provides qualitative data for the ARGs targeted here (presence/absence of a gene in a genus). Therefore, it is important to couple the epicPCR data with 16S rRNA amplicon sequencing data to gain insight into the relative abundance of each detected genus. In addition, several rare genera were detected using epicPCR but not 16S rRNA gene sequencing, which was expected because epicPCR uses different primers to amplify a part of the 16S rRNA gene. In this study, epicPCR provided an overview of the hosts of clinically relevant beta-lactamase genes in Croatian wastewater, which provides a good basis for further research that could be directed towards the cultivation of specific bacterial genera.

Acknowledgement

We thank the staff of WWTP Zagreb and Varazdin for their cooperation in collecting the samples. We acknowledge DNA Sequencing and Genomics Laboratory, Institute of Biotechnology, University of Helsinki for sequencing. We thank CSC—IT Centre for Science, Finland for providing the computational resources.

Author contributions

Svjetlana Dekić Rozman (Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing), Ana Puljko (Investigation), Antti Karkman (Data curation, Methodology, Software, Writing – review & editing), Marko Virta (Conceptualization, Data curation, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing – review & editing), and Nikolina Udiković-Kolić (Conceptualization, Data curation, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing – review & editing)

Conflict of interest

None declared.

Funding

This work was supported by the Croatian Science Foundation (project number IP-2019-04-5539), the Research Council of Finland (Finnish Multidisciplinary Centre of Excellence in Antimicrobial Resistance Research), and European Molecular Biology Organization (Scientific Exchange Grant 9104 awarded to S.D.R.).

Data availability

All 16S rRNA gene and epicPCR sequences have been submitted to the Sequence Read Archive (PRJNA1012856).