Low-level environmental metal pollution is associated with altered gut microbiota of a wild rodent, the bank vole (Myodes glareolus).

Mining and related industries are a major source of metal pollution. In contrast to the well-studied effects of exposure to metals on animal physiology and health, the impacts of environmental metal pollution on the gut microbiota of wild animals are virtually unknown. As the gut microbiota is a key component of host health, it is important to understand whether metal pollution can alter wild animal gut microbiota composition. Using a combination of 16S rRNA amplicon sequencing and quantification of metal levels in kidneys, we assessed whether multi-metal exposure (the sum of normalized levels of fifteen metals) was associated with changes in gut microbiota of wild bank voles (Myodes glareolus) from two locations in Finland. Exposure to increased metal load was associated with higher gut microbiota species diversity (α-diversity) and altered community composition (β-diversity), but not dispersion. Multi-metal exposure and increased levels of several metals (Cd, Hg, Pb and Se) were associated with differences in the abundance of microbial taxa, especially those within the families Clostridiales vadinBB60 group, Desulfovibrionaceae, Lachnospiraceae, Muribaculaceae and Ruminococcaceae. Our data indicate that even low-level metal pollution can affect the diversity of microbiota and be associated with deterministic differences in composition of host gut microbiota in wild animal populations. These findings highlight the need to study a broader range of metals and their cocktails that are more representative of the types of environmental exposure experienced by wild animals.


H I G H L I G H T S
• The effects of metal exposure on wild animal gut microbiota are poorly understood. • We used 16S rRNA amplicon sequencing and measured metal levels in bank vole kidneys. • Multi-metal exposure, and Se were associated with increased gut microbiota α-diversity. • Multi-metal exposure was associated with a shift in microbiota community composition, but not dispersion. • Gut microbiota response to multi-metal exposure and Cd, Hg, Pb and Se is metal and bacterial taxa specific.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Anthropogenic release of metals into the environment is an important cause of stress and detrimental health effects in many wildlife populations (Acevedo-Whitehouse and Duffus, 2009). Biological impacts of exposure to metals are diverse, including effects on metabolism, reproduction and functioning of immune and nervous systems (Sparling, 2017), depending on the metal, dose, organism affected and its developmental timing of exposure. For example, inhabiting an area close to a copper-nickel smelter was associated with a slower growth rate in autumnal moth (Epirrita autumnata) larvae (Eeva et al., 2018), accelerated telomere shortening in great tit (Parus major) but not pied flycatcher (Ficedula hypoleuca) nestlings (Stauffer et al., 2017) and lower levels of antioxidants in Daubenton's bats (Myotis daubentonii) (Ruiz et al., 2019). While diverse effects of exposure to metal pollution on organisms have been quite well-studied, it is less clear whether metal pollution is associated with an impact on the community of host-associated microbes in natural populations.
All animals host communities of microbestheir microbiota. Gut microbiota, the community of microorganisms inhabiting the gastrointestinal tract (Marchesi and Ravel, 2015), can provide essential services to its host such as digestion of otherwise indigestible dietary material (Morrison and Preston, 2016), helping to train the host's immune system (Pickard et al., 2017) and preventing colonization of the gut by potential intestinal pathogens (Buffie and Pamer, 2013). As the community composition of gut microbiota is a key part of host health, it is important to identify determinants of variation in microbiota community composition. For example, host's gut microbiota can be influenced by factors such as host's diet (David et al., 2014;Youngblut et al., 2019), reproductive status (Mallott et al., 2020) and age (Heitlinger et al., 2017).
Changes in gut microbiota composition are also associated with environmental disturbances, for example, land-use change and habitat fragmentation (Ingala et al., 2019), and exposure to environmental pollutants such as radionuclides (Lavrinienko et al., 2018(Lavrinienko et al., , 2020 and pesticides (Yuan et al., 2019). In some cases environmental disturbances have been associated with deterministic changes (i.e., a distinctive shift) in the community composition of microbiota (e.g. in laboratory mice exposed to As (Lu et al., 2014a)). However, a recent study has proposed that certain disturbances can elicit stochastic changes in gut microbiota composition (i.e., an increased dispersion or inter-individual variation) (Zaneveld et al., 2017). This hypothesis has been termed "Anna Karenina principle" (AKP) referring to the first line of Leo Tolstoy's novel Anna Karenina, stating: "Happy families are all alike; every unhappy family is unhappy in its own way". In microbiome research, the AKP predicts that healthy microbiota is characterized by stable microbiota composition with low inter-individual variation, while environmental disturbances cause instability of microbiota, with the gut microbiota of each individual reacting in a different way (Zaneveld et al., 2017).
Laboratory studies have demonstrated that the gut microbiota can be affected by ingestion of various metals, the microbial community's response depending on the type of metal and dose of exposure (Richardson et al., 2018), with host characteristics, such as genotype (Lu et al., 2014b) and sex (Chi et al., 2016) affecting the gut microbiotametal relationship. For example, rats exposed to As, Cr or Ni exhibit a dose dependent decrease in their gut microbiota α-diversity (Richardson et al., 2018), while a decrease in α-diversity was reported in mice exposed to Cu, but not Cd, Pb or Al (Zhai et al., 2017). Exposure to metals can also affect the community composition of laboratory mice gut microbiota, for example in response to As and Cd (Li et al., 2019).
Given the findings of these laboratory studies, it could be hypothesized that changes in gut microbiota would occur in wild animals exposed to metals in their environment. However, most laboratory studies on the potential impacts of metals on host-associated microbiota typically quantify the effects of short-term, single-metal and single-dose exposure (Rosenfeld, 2017). This is in contrast to the experiences of wild animals, that in polluted areas are often exposed to a cocktail of several metals at variable doses (Damek-Poprawa and Sawicka-Kapusta, 2004;Sánchez-Chardi et al., 2009). Moreover, wild animals inhabit complex environments and experience greater pressure to survival (e.g., to find food, care for offspring, avoid predators, and/or defend against pathogens) than their laboratory counterparts. Thus, wild animals might be more sensitive to metals at lower doses, as indicated by the discrepancy in the apparent radiosensitivity of laboratory and wild animals (Garnier-Laplace et al., 2013). Whether exposure to metals at 'environmentally relevant' levels elicits comparable impacts on gut microbiota as observed in laboratory experiments remains to be seen. Despite the importance of gut microbiota to host's health, studies on the effects of low-level metal pollution on wild animal gut microbiota are scarce, and directed towards animals inhabiting polluted aquatic ecosystems, for example, greater flamingos (Phoenicopterus roseus) (Gillingham et al., 2020) and Mongolian toads (Strauchbufo raddei) (Zhang et al., 2016).
Our aim was to quantify how inhabiting an area containing elevated levels of metals affects the gut microbiota of a wild rodent, the bank vole (Myodes glareolus). This widely-distributed species is an excellent model species for studying the effects of pollution on gut microbiota as this animal burrows in soil and has a diverse diet that includes plants, invertebrates, fungi and lichens (Butet and Delettre, 2011), many of which bioaccumulate metals (e.g. see Garty, 2001 for lichens). To address the aim of our study, we tested four hypotheses: (1) exposure to metal pollution is associated with lower α-diversity of gut microbiota; (2) gut microbiota community composition (β-diversity) will differ between animals from low and high polluted areas (3) multi-metal exposure will be associated with higher inter-individual variation in community composition (dispersion, a feature of AKP); and (4) we will be able to identify metals that have notable impact on the abundance of several members of bank vole gut microbiota.
We selected four study areas in the Kemi-Tornio (Fig. 1B) and three study areas in Harjavalta (Fig. 1C). In Harjavalta, the polluted area (H-smelter) was considered to be located within 4 km of the smelter and the unpolluted areas (H-Kokemäki and H-Nakkila) some 8-15 km away from the smelter (Eeva and Penttinen, 2009;Kiikkilä, 2003). In Kemi and Tornio, the expected polluted sites were located less than 8 km from the site of the mine/smelter (Kmine and T-smelter) and the reference areas (K-Simo and T-Kukkola) approximately 17-28 km from the nearest point-polluter (Pöykiö et al., 2005a(Pöykiö et al., , 2005b. We set 3-10 trapping sites within each study area ( Fig. 1B-C) in which 12 Ugglan Special No. 2 live traps (Grahnab, Sweden) baited with sunflower seeds and potato were placed, in three grids that were at least 15 m from the nearest road (4 traps per grid, an inter-trap distance of 10-15 m, and an inter-grid distance of 50 m). Traps were set in the late afternoon and checked the following morning for two consecutive nights. Adult voles were euthanized by cervical dislocation and immediately transferred to dry ice and stored at −80°C until further processing. Bank vole (n = 230) head width (a proxy for age, Kallio et al., 2014), weight, sex and gravidity status were determined during dissection. Organs and faecal samples from the distal 2 cm of the colon were immediately put on dry ice and stored at −80°C. All animal procedures were conducted in accordance with the Finnish Act on the Use of Animals for Experimental Purposes, with permission from the Finnish Animal Experiment Board (ESAVI-3981-2018).

Determining concentration of metals in bank vole tissue
We selected a subset of 84 animals (41 female, 43 male) from each study area (12 animals per study area) to determine levels of 71 elements in the right kidney of bank voles.
Chemical analyses were made by ALS Scandinavia (Luleå, Sweden, https://www.alsglobal.se/en) by inductively coupled plasma-sector field mass spectrometry (ICP-SFMS, ELEMENT2, ThermoScientific, Bremen, Germany) using a combination of internal standardization and external calibration. Samples were prepared using closed vessel, microwave assisted acid digestion. Handling of all samples and digests was done in clean laboratory areas (Class 10,000), following necessary precautions to reduce contamination risks (Rodushkin et al., 2010). The ERM-BB186 Pig Kidney certified reference material (JRC, Belgium) was prepared and analysed together with vole samples and used for assessment of accuracy and precision of the analytical method. Further details on preparation method, ICP-SFMS instrumental parameters, operation conditions and figures of merit can be found elsewhere (Engström et al., 2004). An average element concentration was used for samples with replicate measurements available, and the elements with levels below the limit of detection were transformed to half of the corresponding limit. Metal levels (μg/g dry weight) for all animals can be found in Table S1. To study the effects of metal pollution on gut microbiota we selected a subset of 13 metals (Cd, Co, Cr, Cu, Fe, Hg, Mn, Mo, Ni, Os, Pb, V, Zn) and two metalloids (As, Se), as these elements either have known impacts on wildlife in the study areas (e.g. Rodushkin et al., 2011;Stauffer et al., 2017) or have apparent impacts on rodent health and host gut microbiota (Duan et al., 2020). For brevity, we refer to these 15 elements collectively as "metals" throughout the manuscript.

Faecal DNA extraction and 16S rRNA amplicon sequencing
Total DNA was extracted from bank vole faeces using QIAGEN PowerFecal DNA Kit (QIAGEN, Germany) following the manufacturer's protocol (except step 13using 600 μl of supernatant). Two hundred samples (i.e., excluding samples with low DNA yields) were sent to Institute for Molecular Medicine Finland (FIMM, https://www.helsinki. fi/en/hilife-fimm) for amplicon sequencing of the V4 region of the 16S ribosomal RNA (rRNA) locus using the 515F/806R primers (Caporaso et al., 2011) on an Illumina MiSeq (250 bp paired-end reads). Bank vole metadata used in the analyses are provided in Table S2.

Sequencing data processing
Reads were demultiplexed at the sequencing facility (FIMM). We used cutadapt v.1.10 (Martin, 2011) to remove Illumina adapter sequences. Demultiplexed, paired-end reads (n = 17,518,994) were processed using QIIME2 v.2019.10 (Bolyen et al., 2019), with the DADA2 plugin (Callahan et al., 2016) used to remove primers, truncate (reverse reads at 176 bp due to a quality drop) and merge paired-end reads, remove chimeric sequences and to cluster the sequences into amplicon sequence variants (ASVs) (Callahan et al., 2017). This procedure generated 10,592,451 sequences and 5287 ASVs in 193 samples (seven samples failed during library preparation and sequencing), with an average of 54,192 reads and 233 ASVs per sample. We filtered out low abundance (<10 reads in the entire data set) ASVs, and assigned taxonomy to the representative sequences using a pre-fitted sklearn-based classifier trained on the SILVA 132 database from 515F/806R region of the 16S rRNA gene (Yilmaz et al., 2014) using a QIIME2 plugin feature-classifier (Bokulich et al., 2018;Pedregosa et al., 2011). After removing ASVs (n = 26) that could not be assigned to a bacterial phylum, and the ASVs classified as mitochondria or chloroplasts, we created a mid-point rooted phylogenetic tree using FastTree (Price et al., 2010). This final dataset had a total of 10,582,131 reads and 3653 ASVs and an average of 54,119 (range 9704-85,330) reads and 224 ASVs per sample. Samples were rarefied to an even depth of 27,000 reads per sample (leaving n = 192 samples), and these rarefied data were used for analyses unless otherwise stated. Raw reads are publicly archived in NCBI Sequence Read Archive under BioProject number PRJNA702897.

Statistical analyses
Statistical analyses were completed in R v.4.0.2 (R Core Team, 2020) except when stated otherwise. The relative physiological condition of animals was estimated using a body condition index (BCI), calculated as the standardized residuals of a linear regression of body mass against head width, where positive BCI values indicate better condition and possibly higher energy reserves (Koskela et al., 2004;Schulte-Hostedde et al., 2001). Animals that were gravid (n = 18) or did not have complete metadata (n = 3) were excluded from all analyses.

Metal levels and pollution in study areas
To determine the metal burden of animals, we calculated the total metal load (TML) of each animal as the sum of normalized (x -min(x) / (max(x) -min(x))) levels (Santamaría et al., 2012) of the 15 metals. Normalization of metal level was necessary to account for variation in the scale of metal levels. To assign animals to a pollution group (high or low) we fitted linear model with TML of individuals (n = 83, excluding one gravid animal) as the response variable and study area, bank vole head width, BCI and sex as the explanatory variables. Based on the results of linear model (Table 1) there were significant differences in TML between study areas -animals from study areas H-Kokemäki, K-Simo and K-mine had significantly lower TML than animals from study areas H-Nakkila, H-smelter, T-smelter and T-Kukkola. Animals from the same study area were assigned to the same pollution group (Table 1), based on limited (usually <1 km during breeding season) dispersal of bank voles (Kozakiewicz et al., 2007). High pollution level was characterized by significantly (p < 0.05) higher levels of As, Co, Cr, Fe, Mo, Ni, Os, Pb, Se, and V, for example the mean level of Cr was 1.4 and As 5.2 times higher in animals from study areas in group the "high" compared to "low" (Fig. 2). Metal levels between all study areas were compared by Kruskal-Wallis tests using Dunn's test of multiple comparisons with Benjamini-Hochberg correction (Table S3) via package dunn.test v.1.3.5.

Gut microbiota analysis
Effects of metal pollution on gut microbiota of bank voles was examined (1) among the two pollution groups, to examine the possible effects of multi-metal exposure in 173 animals and (2) in the subset of 80 animals for which we had data on individual metal concentrations. Association of BCI and multi-metal exposure and individual metals was examined in 179 and 83 animals, respectively (including individuals that failed sequencing). Alpha diversity was calculated as (1) Shannon's diversity index and (2) observed ASV count using phyloseq v.1.32.0 (McMurdie and Holmes, 2013), and as (3) Faith's phylogenetic diversity using picante v.1.8.2 (Kembel et al., 2010). Beta diversity was assessed using Bray-Curtis dissimilarity and UniFrac (Lozupone and Knight, 2005) distance which were calculated using phyloseq. We used the Songbird v.1.0.1  plugin in QIIME2 to identify ASVs most associated with multi-metal exposure and individual metals. Songbird uses multinomial regression models on unrarefied datasets to perform differential ranking analysis. This analysis produces differentials ("logarithm of the fold change in abundance of a taxa between two conditions" ) which can be ranked based on their association with the covariate of interest. For all differential abundance analyses we selected the top 10% ASVs as numerators and bottom 10% as denominators (reference frame) since these ASVs have the strongest positive or negative associations with variables of interest (pollution group and individual metals) . The resulting differentials were visualized and log-ratios of top and bottom ASVs were computed using Qurro v.0.7.1 (Fedarko et al., 2020) and have been deposited in Figshare (Brila, 2021a).
2.5.2.1. Association of multi-metal exposure with gut microbiota. To assess the association of multi-metal exposure with gut microbiota α-diversity, we fitted linear mixed models (LMMs) with each α-diversity metric as response variable and pollution group as the explanatory variable, and included BCI, sex, head width and location as covariates in all models to account for the effects of host trait and environmental heterogeneity. To evaluate if multi-metal exposure is associated with distinct shifts in community composition we used permutational multivariate analysis of variance (PERMANOVA) with 999 permutations with sequential evaluation of terms via adonis in vegan v.2.5-6 (Anderson, 2001). Each β-diversity metric was used as a response variable, and pollution group as the explanatory variable, with host traits (head width, BCI and sex) and environmental traits (location and study area) included as covariates in the same model. The level of microbiota community dispersion (inter-individual variation) between pollution groups was examined using vegan functions betadisper and permutest. Core microbiota, defined as ASVs present in >50% of samples, was determined for each pollution group using microbiome v.1.10.0 (Lahti and Shetty, 2020). Association of multi-metal exposure with BCI was assessed using LMM with the pollution group as explanatory variables and host traits (head width and age) as covariates. ASVs associated with multimetal exposure were found using optimized multinomial regression model (differential prior 0.3, epochs 10,000, random test samples 17) with pollution group as the explanatory variable (low pollution group used as the reference) and host traits (head width, BCI, sex) and location included as covariates. The multinomial regression model uses only ASVs present in a minimum of 10 samples (resulting in 751 differentially ranked ASVs) and was compared to the baseline (covariates only) model (Q2 = 0.06) and null model (Q2 = 0.12). The log-ratios were calculated using top 10% ASVs (n = 75) as numerators and bottom 10% (n = 75) as denominators (reference frame). To statistically verify these findings, the log-ratios were used as the response variable in LMM (n = 172, 1 sample had log-ratio containing zero), with location and host traits included as covariates.
2.5.2.2. Association of individual metals with gut microbiota. Effects of individual metals on α-diversity metrics and BCI were quantified using LMMs with levels of 15 metals as explanatory variables while controlling for host traits (sex and head width, and for α-diversity, BCI) and location. Prior to fitting LMMs, metal levels were log 10 transformed to minimize variation in skewness and scaling. Polynomial terms were included based on non-linear relationships detected using exploratory plotting. Best set of models were identified based on AICc (Akaike Information Criterion corrected for small sample sizes) (Akaike, 1974;Sugiura, 1978) calculated using the pdredge function in MuMIn v.1.43.15 (Bartoń, 2020) via R v.3.6.3. We calculated full (unconditional) averaged estimates from all models within 2 ΔAICc from the model with the lowest AICc (Burnham et al., 2011;Symonds and Moussalli, 2011) to estimate the effects of metals on BCI and α-diversity metrics using MuMIn v.1.43.17. To identify which metals were most associated with changes in gut microbiota composition, we used bioenv in vegan v.2.5-6 (Oksanen et al., 2019) to find the metals with highest Spearman's rank correlation with β-diversity metrics. We then performed partial-Mantels tests between identified metals (Hg and Os) and Bray-Curtis and UniFrac distances based on Spearman's rank correlation with 9999 permutations using mantel.partial function in vegan, while controlling for the effects of host traits, location and study areas. Multinomial models were used to identify ASVs associated with specific metals, with host traits (head width, BCI, sex), location and study area included as covariates in all models. Final optimized (differential prior 0.5, epochs 10,000, random test samples 8) model included Cd, Hg, Pb and Se, and was compared to the baseline model (Q 2 = 0.08) and null model (Q 2 = 0.17). The final model produced differential ranking of 435 ASVs. Log-ratios for individual metals were calculated as described previously, selecting top 10% and bottom 10% of ASVs (n = 43) to calculate the log-ratios. To test these findings statistically, we used log-ratios as the response variable in LMMs (four separate models), with location and host traits included as covariates.

Pollution levels in study areas
Metal levels quantified in bank vole kidneys were overall higher than the concentrations of metals reported in bank voles collected from five long-term monitoring sites in Sweden, where the distance to nearest point polluter in all but one sites is >10 km (Ecke et al., 2020) and comparable or lower than those found in rodents from other polluted sites in Europe (Damek-Poprawa and Sawicka-Kapusta, 2004;Sánchez-Chardi et al., 2009). High metal load at T-Kukkola, which is located >17 km from the nearest obvious source of pollution (Fig. 1B), is likely due to air transport of fine particles from the ferrochrome production plant (Bulut et al., 2009;Rodushkin et al., 2007). The high metal load at H-Nakkila might be a consequence of the study area being downstream the river Kokemäenjoki that runs through Harjavalta, where several tons of nickel and 1.3 t of cobalt were spilled during an accident in July 2014 (KVVY ry, 2016). Indeed, we found elevated levels of Ni and Co in the kidneys of bank voles from study areas H-Nakkila and H-smelter (Fig. S1). For example, the mean level of Ni was 40% and the mean level of Co 37% higher in H-Nakkila and H-smelter study areas, respectively, when compared to the upstream H-Kokemäki study area. Relatively low levels of metals detected in kidneys of bank voles at the Kemi mining area (K-mine) are likely a consequence of the mine at Kemi being underground (since 2006), and that most of the initial processing of ore also occurs underground or in closed facilities. Hence, the amount of dusting of the surrounding environment is minimised (Kauppila et al., 2011). The ore (chromite or chromium oxide) mainly contains a trivalent form of Cr that is less biologically available, and thus is considered non-toxic (unlike hexavalent Cr) to humans and wildlife (Beukes et al., 2017).

Association of multi-metal exposure and gut microbiota diversity and composition
Inhabiting areas with elevated levels of metals was associated with higher α-diversity and some differences in β-diversity of gut microbiota. Contrary to our hypothesis, all metrics of α-diversity tended to be higher in the gut microbiota of animals inhabiting areas belonging to the high pollution group (Fig. 3A-C), although these differences were statistically significant only for Faith's phylogenetic diversity index (p = 0.041) (Table S6).
We found significant differences in community composition (β-diversity) between high and low pollution groups (PERMANOVA, Bray Curtis dissimilarity R 2 = 0.02, p = 0.001 and UniFrac distance metric R 2 = 0.01, p = 0.001) (Fig. 3D-E). Additionally, community composition was affected by host traits (BCI and sex), and environmental factors (location and study area) ( Table S7). The observed differences in community composition between locations however, could be a consequence of significantly different dispersion in both Bray-Curtis dissimilarity (permutest pseudo-F = 12.031, p = 0.002) and UniFrac metric (permutest pseudo-F = 6.194, p = 0.01) as it may influence the results of PERMANOVA (Anderson and Walsh, 2013). Counter to our hypothesis, metal pollution was not associated with greater community dispersion (Bray-Curtis permutest pseudo-F = 3.475, p = 0.071, UniFrac permutest pseudo-F = 0.331, p = 0.559). There were no significant (p > 0.05) differences in BCI between pollution groups (Table S8).

Individual metal associations with gut microbiota diversity and composition
Few metals were associated with altered gut microbiota α-diversity or community composition. While no metal had a significant association with Shannon diversity and Faith's phylogenetic diversity, an increase in Se level was positively associated with (model average relative variable importance = 1) the number of observed ASVs (Table S9). An increase in Hg and Os concentrations was associated with altered community composition of gut microbiota (Bray-Curtis dissimilarity r = 0.21, p = 0.002 and r = 0.26, p = 0.001 and UniFrac distance metric r = 0.26, p ≤0.001 and r = 0.18, p = 0.008, respectively). Neither Hg nor Os were associated with changes in bank vole body condition. However, high levels of As and Cr were negatively associated with the BCI of bank voles, while Co and Zn showed a positive association, and Mo a non-linear relationship with BCI showing a dose dependent positive, hump-shaped relationship (Table S10).

Association of differentially abundant ASVs with pollution and specific metals
We found significant differences between the top and bottom ranked ASVs associated with multi-metal exposure (p = 0.003, Table S11). Some of the ASVs positively associated with multi-metal exposure include members of genera Anaerostipes, Odoribacter and Lachnospiraceae NK4A136 group, while ASVs negatively associated with multi-metal exposure include ASVs from genera Desulfovibrio, Treponema 2 and Ruminococcus 1 (Table S12). However, we found that many families had ASVs that were positively and ASVs negatively associated with multi-metal exposure. For example, Clostridiales vadinBB60 group had 9 ASVs positively and 3 ASVs negatively associated with multi-metal exposure, while Lachnospiraceae had 20 positively and 22 negatively associated ASVs. Furthermore, variation in response to multi-metal exposure was evident at the genus level, for example within family Ruminococcaceae ASVs from genus Ruminococcaceae NK4A214 group were negatively associated with multi-metal exposure. However, genera Ruminiclostridium 6 and Ruminiclostridium 9 had ASVs that were positively and ASVs that were negatively associated with multi-metal exposure. Five bacterial families: Lachnospiraceae (n = 42), Muribaculaceae (n = 24), Desulfovibrionaceae (n = 22), Ruminococcaceae (n = 15) and Clostridiales vadinBB60 group (n = 12) contained 76.7% of all ASVs identified as either positively (top 10% differentially ranked ASVs) or negatively (bottom 10% differentially ranked ASVs) associated with multi-metal exposure.
For Cd, Hg, Se and Pb (p < 0.001 for all models, see Table S11), we identified metal specific associations with different ASVs. For example, Hg was negatively but Se positively associated with ASVs from genera Ruminococcus 1 and Ruminiclostridium 5, while only Se was positively associated with Intestinimonas ASV (Table S13). Many of the differentially abundant ASVs identified as top or bottom-ranked for individual metals (Cd, Hg, Se or Pb) were from the same genera, for example, Alistipes, Desulfovibrio and Ruminiclostridium. Furthermore, some ASVs showed similar response to multiple metals, for example, increased levels of Cd, Hg and Pb were negatively associated with Anaerostipes ASV. Five bacterial families were associated with changes in levels of all four metals -Clostridiales vadinbb60 group, Desulfovibrionaceae, Lachnospiraceae, Muribaculaceae and Ruminococcaceae (Table S13). Though all five families had ASVs positively and ASVs negatively associated with individual metals, we found that some genera within these families had strong positive or negative associations. For example, within family Ruminococcaceae, ASVs from genera Oscillibacter and Ruminiclostridium 9 were negatively associated with increased levels of Cd, while increased levels of Se were positively associated with Oscillibacter and Ruminococcus 1 ASVs. However, even within genera we saw contrasting response of ASVs, for example, Lachnospiraceae NK4A136 group had ASVs ranked as most and as least associated with increased levels of Hg and Pb, indicating that the response to individual metals is genus and even ASV specific.

Discussion
Gut microbiota is thought to be an important part of wildlife health, but little is known about the association between anthropogenic environmental disturbance and gut microbiota diversity and composition. Here, we found that exposure to low-level chronic metal pollution was associated with altered gut microbiota of wild bank voles. Animals inhabiting areas with elevated levels of metals had higher species diversity (α-diversity), but there was no change in the level of community dispersion. We found a subtle shift in community composition (β-diversity), while bacterial families most associated with multi-metal exposure included Clostridiales vadinbb60 group, Desulfovibrionaceae, Lachnospiraceae, Muribaculaceae and Ruminococcaceae. Metals most associated with altered abundance of several gut bacteria taxa were Cd, Hg, Pb and Se. Hence, even low-level metal pollution associates with altered gut microbiota of wild mammals.

Association of metal pollution and alpha diversity of gut microbiota
Our study revealed a trend of higher α-diversity and a more diverse core microbiota in the high pollution group animals. However, only Faith's phylogenetic diversity was statistically different between the pollution groups, with higher divergence of phylogenetic lineages seen in high pollution group animals. This trend was unexpected as studies on laboratory rats and mice have found metal exposure to be associated with either no differences in or reduced α-diversity of gut microbiota (Richardson et al., 2018;Zhai et al., 2017). However, some studies on humans inhabiting areas with metal pollution have shown an comparable increase in α-diversity of their gut microbiota (Eggers et al., 2019;Shao and Zhu, 2020).
Variation in microbiota α-diversity in response to exposure to metal pollution in laboratory and field studies might be an example of the intermediate disturbance hypothesis (IDH, Connell, 1978) in gut microbiota. IDH predicts that species diversity is expected to be highest at intermediate levels of disturbance, explained by competitioncolonization trade-off wherein less competitive species are able to reproduce and colonize the environment before being eliminated by more competitive species, resulting in a coexistence of the two with higher diversity of the community (Sheil and Burslem, 2013). Thus, it might be plausible that low-intermediate metal exposure allows coexistence of many species of gut microbiota, whereas no pollution would result less diverse and more stable, "core" microbial community, while high levels of metal pollution would favour metal-resistant microbial taxa. While the IDH has been examined in freshwater microbial microcosms (Gibbons et al., 2016), it is not I. Brila, A. Lavrinienko, E. Tukalenko et al. Science of the Total Environment 790 (2021) 148224 known whether IDH applies to gut microbiota of animals exposed to environmental pollutants. Further experiments that quantify biological impacts (single metal and multi-metal) across a wide pollution gradient are needed to understand the applicability of IDH to different types of disturbances and microbiota. The positive association between α-diversity and Se is intriguing, given the apparent increase in gut microbiota α-diversity following Se supplementation in laboratory rodents (Cheng et al., 2020;Kasaikina et al., 2011). There is some indication that Se is important for gut microbiota and host health, as this essential element is involved in the maintenance of gut barrier function and resistance to pathogen infection (Zhai et al., 2018). Though essential, Se can be toxic at high levels. The levels of Se in kidneys of bank voles observed in our study (mean 2.9, range 1.68-5.23 μg/g, dry weight) were lower than those found in the liver of three rodent species inhabiting an area with elevated levels of Se and comparable to levels found in reference areas (Santolo, 2009), but higher than those found in kidneys of wild bank voles in Sweden (mean: 0.69, range 0.18-1.99 μg/g, dry weight; Ecke et al., 2020). However, as animal species vary in their sensitivity to metals (Allard et al., 2010) whether the association between Se and α-diversity is beneficial or detrimental to the host remains unclear.
The significance of higher gut microbiota α-diversity in areas with greater total metal load is unclear. On the one hand, α-diversity has been associated with greater functional potential of bacterial communities, and resistance to environmental perturbations and pathogen invasion (Lozupone et al., 2012). On the other hand, high α-diversity might be disadvantageous (Reese and Dunn, 2018), for example, with the microbial community becoming less stable and potentially containing more opportunistic bacteria (Coyte et al., 2015). The host might, therefore, benefit from gut microbiota α-diversity staying within a certain range. What such range may be for wild bank voles and whether the observed differences in α-diversity are sufficient to elicit, for example, a significant change in microbiota functions or if the differences represent a change in the number of opportunistic pathogens warrants future investigation.

Low-level metal pollution and gut microbiota community composition and dispersion
Though low-level metal pollution was associated with a slight but significant shift in community composition, the increase in level of community dispersion, that is one feature of the "Anna Karenina principle" (Zaneveld et al., 2017), was not observed. While the variance in community composition explained by metal pollution was low, the effect of metal exposure was comparable to the effect of location, and greater than that of host traits (BCI and sex). These findings indicate that lowlevel metal pollution is associated with slight deterministic rather than stochastic changes in gut microbiota community, with individuals showing a similar, and potentially predictable response to multi-metal exposure. The high inter-individual variation in community composition seen in our study is a feature of microbiota studies, especially in wild animal gut microbiota (Viney, 2019) that have diverse backgrounds and developmental experiences. Thus, we expect that in addition to metal pollution, bank vole microbiota will be affected, for example, by variation in host genotype (Suzuki et al., 2019) and diet (Youngblut et al., 2019), that are features that are not readily measured or controlled in studies of wild populations.
While higher microbiota community dispersion is a hypothesized indicator of lower community stability, few studies on the effects of environmental pollution on gut microbiota have explicitly tested for differences in community dispersion. Nevertheless, a shift in community composition and lower dispersion is apparent in the gut microbiota of laboratory mice exposed to Cd but not As (Li et al., 2019) and in wild bank voles exposed to environmental radionuclides in Chernobyl (Lavrinienko et al., 2020) and an increase in dispersion was associated with long-term multi-metal exposure in humans inhabiting mining and smelting areas (Shao and Zhu, 2020). Thus, the type of gut microbiota community response (stochastic or deterministic), if any, depends on the type, duration and magnitude of pollution and is not a clear marker per se of environmental disturbance.
Though Hg and Os did not associate with a change in the body condition index (a proxy for relative physiological condition) of bank voles, the observed association of these metals with the community composition of bank vole gut microbiota is relevant. Potential effects of Hg on the gut microbiota are far less studied than those of As, Cd and Pb (Duan et al., 2020). Furthermore, no study so far has addressed whether gut microbiota may be altered when animals are exposed to Os -a highly toxic metal, impacts of which are rarely studied in wildlife, but see (Rodushkin et al., 2011) and (Ecke et al., 2018) for potential somatic effects in bank voles. Conversely, the negative association of BCI with As and Cr, and dose dependent association between BCI and level of Mo, indicates that low-level exposure to these metals might impact the physiological condition of wild bank voles, yet these metals were not associated with changes in community composition. Similarly, exposure to Cr in laboratory rats was associated with decreased weight of kidneys and liver, but had little effect on community composition of gut microbiota (Richardson et al., 2018). This highlights that in multi-metal exposure individual metals might be able to affect different components of host health, e.g., affecting host directly (e.g., decreased body condition index or organ weight) or indirectly by affecting host's gut microbiota.

Metal pollution is associated with varying response of bacterial taxa
Composition of bank vole gut microbiota being dominated by members of the Bacteroidetes, Firmicutes and Proteobacteria is similar to the taxonomic composition reported in bank voles in other studies (Lavrinienko et al., 2018(Lavrinienko et al., , 2020. Three of the five families, namely Ruminococcaceae, Lachnospiraceae and Muribaculaceae that were most associated with metal pollution are common and abundant members of bank vole gut microbiota (Lavrinienko et al., 2018(Lavrinienko et al., , 2020. Thus, it might not be surprising that these families had the most ASVs positively or negatively associated with multi-metal exposure and individual metals. Consistent with our findings of varying response within bacterial families and genera, previous studies on metal pollution and gut microbiota have reported contrasting findings on the effects of metal pollution on several bacterial families. For example, most studies report negative association of metals with Lachnospiraceae (Duan et al., 2020), however an increase was reported in the gut microbiota of humans exposed to multiple metals in the vicinity of a Pb\ \Zn mine (Shao and Zhu, 2020). We were able to identify some bacterial genera with uniform, metal-specific response. For example, the negative association of several Ruminococcaceae genera with Cd is comparable to findings of other studies, e.g., decrease in relative abundance of Ruminococcaceae observed in mice exposed to Cd and Pb (Zhai et al., 2017). Furthermore, we detected a positive association of multiple Ruminococcaceae genera with increased levels of Se, an element that has been linked to mitigation of toxicity of several metals including As and Cd (Zwolak, 2020).
The complexity of the gut microbiota response is apparent, given that the response to multi-metal exposure and individual metals was genus and ASV specific. The observed patterns of differentially abundant taxa suggest that using a specific family of bacteria as a marker of environmental disturbance of the gut microbiota is likely to be misleading. Furthermore, whether the observed changes in abundance of bacteria affect the services provided to the host remains unclear, as functional redundancy of gut microbiota has been reported, with same services provided by multiple bacterial taxa (Heintz-Buschart and Wilmes, 2018).
In this study, we demonstrate that exposure to low-level metal pollution is associated with altered gut microbiota of wild bank voles. The trend of higher α-diversity associated with multi-metal exposure shows similarities with environmental exposure to multiple metals in humans and differs from the findings of single metal exposure in laboratory animals. Thus, it seems likely that exposure to multiple metals simultaneously is associated with a different microbiota response than single-metal exposure and suggests that the effects of environmental metal exposure on gut microbiota of wild animals cannot be extrapolated from laboratory studies. Though we identified several bacterial families associated with multi-metal exposure and with elevated levels of Cd, Hg, Pb and Se, the observed variation in response within and among families highlights the complexity and diversity of effects of metal exposure on gut microbiota. Future work, potentially incorporating analysis of host diet and health, should aim to examine the functional significance and potential health effects of these subtle changes in gut microbiota diversity and composition.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.