Perch was sampled at two locations (Kalmar, 56°40.306′N; 16°21.578′E and Figeholm, 57°22.321′N; 16°33.340′E) separated by a swimming distance of approximately 80 km along the southeast Baltic coast of Sweden. This spatial separation exceeds the general dispersal pattern and home ranges (∼20 km) of Baltic Sea perch (Craig, 2000; Ahlbeck Bergendahl et al., 2017; Hansson et al., 2019) such that the locations likely harbor specific host populations with non-overlapping home ranges. This is also supported by comparisons of microsatellite data from Baltic Sea perch, including our study area, that show genetic clusters (populations) at a much finer spatial resolution than our two locations (Bergek and Björklund, 2009; Olsson et al., 2011). The spatial separation of populations also implies exposure to different environmental conditions that may shape the microbial community, such as sediment load, bottom substrates, and habitat heterogeneity (Caporaso et al., 2011; Walters and Martiny, 2020).

The Figeholm population resides in a habitat that consists of an approximately 5 km wide archipelago with distinct depth gradients as well as two streams discharging into the area. The drainage area is sparsely populated (<1,000 inhabitants), semirural and covered by coniferous forest suggesting limited direct anthropogenic impacts on the aquatic habitat. The Kalmar perch population inhabit a coastal area that generally lacks adjacent archipelago and freshwater inflows. The sample location is in the heart of the urban area of Kalmar municipality (>50,000 inhabitants) with substantial man-made modifications of the habitat including embankment, beachfront buildings and extensive hard surfaces resulting in direct anthropogenic impacts on the habitat and water quality.

Focal individuals [Kalmar: n = 30; Figeholm: n = 9; average size of 33.3 ± 3.6 cm (mean ± SD)] were captured between the 9th and 16th of October 2013. To avoid cross-contamination among individual fish microbiomes, individuals were captured one-by-one using rod and reel fishing. Microbiome samples were obtained from fish immediately after capture following rinsing with sterile MQ-water to eliminate bacterial cells associated with the water column. Pre-defined spots (approximately 2 × 2 cm) on the fish body were sampled with a sterile cotton swab that was twirled four times on the spot and then transferred to an Eppendorf-tube with 750 μl TE buffer (Tris-EDTA, 10:1). To minimize cross-contamination among samples all equipment used for sampling was sterilized with 70% ethanol between each sample. All samples were kept on ice until being stored at −80°C. The majority (32 of 39) of the fishes were females. Immediately after sampling, each fish was stunned and terminated by cervical dislocation. All applicable national guidelines for the care and use of animals were followed. Ethical approval for the study was granted by the Ethical Committee on Animal Research in Linköping, Sweden (Dnr. 33–14 and 10–14).

We obtained two samples (dorsal and ventral) from the right side of every individual (n = 39). From a subset (n = 16) of these individuals, representing both populations, we also took two additional corresponding samples from the left side of the body (both dorsal and ventral) to be able to evaluate heterogeneity and sample repeatability within individuals (Figure 1). Data thus consists of 110 samples that represent 2–4 samples from 39 fish individuals originating from 2 different populations. Specifically, 82 samples from the Kalmar population of 30 fish individuals that were caught on three occasions (9, 11, and 15 October 2013; 4 samples from 11 individuals and 2 samples from 19 individuals) and 28 samples from the Figeholm population distributed among 9 fish individuals that were caught on one occasion (16 October 2013; 4 samples from 5 individuals and 2 samples from 4 individuals).

At each sampling occasion and location, surface water samples (1 L) were taken to enable comparisons of the microbial communities in the water column with those present on the fish skin. Water samples (n = 4) were vacuum filtered through a Supor membrane filter (Pall Corporation, pore size 0.22 μm, Ø 47 mm). The filter was transferred to a sterile tube containing 1.8 mL TE-buffer (Tris-EDTA, 10:1) and stored at −80°C until DNA extraction.

DNA was extracted from the microbiome samples collected with sterile swabs (n = 110) using the QIAmp DNA Mini Kit (QIAGEN, Germany) starting from step 3 in the protocol by the manufacturer. Final elution volume was set to 100 μl to increase the DNA yield. To minimize cross-contamination of samples, the nipper used to take out swabs was sterilized with 70% ethanol between each sample. DNA from water filters (n = 4) was extracted with the DNeasy PowerWater kit (Qiagen, Germany) and obtained concentrations were measured using NanoDrop 2000. Sequencing libraries were prepared by using the 16S rRNA gene primer pair 341F and 805R (Herlemann et al., 2011) following the PCR-protocol by Hugerth et al. (2014) modified by Lindh et al. (2015), with the exception of an additional five cycles in the first PCR. The process of adding Illumina adapters and index sequences was conducted according to Lindh et al. (2015). Final DNA concentrations were measured with a Qubit^® 2.0 Fluorometer. The resulting purified (individually barcoded) amplicons were pooled into three libraries that were purified with gel purification kit (E.Z.N.A.^® Gel purification kit, Omega Bio-Tek, Inc.) and sequenced on three separate occasions (n = 70 and 40 on each occasion for microbiome samples, and n = 4 for water samples) on the Illumina MiSeq platform (Illumina, United States) with 2 × 300 bp paired-end settings at Science for Life Laboratory (SciLifeLab, Stockholm, Sweden).

Both source population and individual accounted for variation in the microbiome composition, but host individual explained a larger proportion of the total variation than did population (RDA, n = 110, effect of population: F_{1, 108} = 1.80, P = 0.007, R² = 0.007; effect of individual: F_{38, 71} = 1.76, P = 0.001, R² = 0.21; Figures 5A,B and Supplementary Figure 3). The amount of variation in microbial community composition was comparable between the two populations (ANOVA, effect of population on Euclidean distance to centroid: F_{1, 108} = 0.046, P = 0.83), indicating that the difference in microbiome composition was not due to a difference in variance among samples within populations. When analyzing host populations separately, individuals within populations were still significantly different from each other (RDA, effect of individual, Figeholm: F_{8, 19} = 2.10, P < 0.001, R² = 0.25; Kalmar: F_{29, 52} = 1.65, P < 0.001, R² = 0.19).

The total number of phyla discovered among all fish skin microbiomes was 57, of which 52 were present in the Kalmar and 47 in the Figeholm population, respectively. Comparisons of the 12 most abundant microbiome phyla in each host population showed that 8 of these 12 phyla were shared between the two populations, but their prevalence varied greatly among individuals (Figures 5A,B). Moreover, these eight phyla included three phyla that were not detected in the water samples: Euryarchaeota, Deinococcus-Thermus, and Gemmatimonadetes. A few phyla had notably high mean relative abundance in one host population while being rare in the other, and this was attributable to variation among individuals. For example, the phyla Tenericutes (majority of sequences belonging to family Mycoplasmataceae) had mean relative abundances 28.7% ranging between 0.007 and 88% in Figeholm, but 4.4%, varying between 0.001–35% in Kalmar, and Patescibacteria (majority of sequences belonging to family Gracilibacteria bacterium JGI 0000069-P22) had mean relative abundance 14.4%, with a range of 0.002–80% in Figeholm while the mean relative abundance was 0.8%, ranging between 0.001 and 9.2% in Kalmar (Figures 5A,B). Within population comparisons revealed that except for the phyla present in all samples from both populations (Figure 4), Patescibacteria was present in all samples from the Kalmar host population. No other phyla were present in all samples from either of the two populations, respectively.

We identified a total of 1,278 ASVs in the bacterioplankton communities from water samples, distributed among 16 and 12 phyla in Kalmar and Figeholm, respectively. None of these phyla were exclusively found in water. The most abundant taxonomic groups, Actinobacteria, Bacteroidetes, and Proteobacteria, were shared between water samples from the two geographic locations, [overall mean abundance 29.4% (range = 18.4–58.7%); 9.9% (range = 4.3–21.1%); 9.1% (range = 4.8–18.2%)] (Figure 5C). These phyla were also the most abundant in skin microbiomes but included more clades at lower taxonomic levels (23, 13, and 87 vs. 13, 8, and 31 orders in microbiome and water samples respectively Supplementary Data File 2). The community composition of bacterioplankton communities in the water was significantly different from that in the skin microbiomes (RDA, effect of sample type: F_{1, 112} = 17.4, P < 0.001, R² = 0.13; Figure 5).

Species richness was highly correlated between left and right samples for both dorsal and ventral body sites within individuals (dorsal: 85.6%, F_{15, 16} = 12.9, P < 0.001; ventral: 87.1%, F_{15, 16} = 14.5, P < 0.001, Figure 6). The microbiome community composition did not differ significantly between samples taken from the dorsal and ventral body parts of the same host individuals, but the significant variation among hosts as reported above (see the effects of individual in the RDA analyses in the previous section) was further supported (PERMANOVA random effect of host individual: F_{15, 47} = 2.16, P < 0.001, R² = 0.22; effect of body site within host: F_{1, 47} = 1.11, P = 0.30, R² = 0.03, Figure 6). There was a significant difference in variance among samples taken from the same host individual (PERMDISP: F_{15, 48} = 10.2, P < 0.001, Figure 2). This reflected that the microbiome community composition was more heterogeneous (i.e., differed more between body sites) in some host individuals than in others (note individual 25 and 26 in Figure 2).

We found more ASVs at the population than individual level (Figure 3), indicating that microbiome composition varied among individual fish hosts. When examining the prevalence of phyla, only three were present in all samples: Actinobacteria, Bacteroidetes, and Proteobacteria, and their relative abundances were highly variable among samples (Figure 4). These phyla were also present in the bacterioplankton communities in the water. However, their relative abundances were much higher in the water compared to the microbiome samples (Figure 5), and the composition at lower taxonomic levels was far more diverse in fish skin microbiomes (Supplementary Data File 1). Other studies have also found the three aforementioned taxa to be the most prevalent in fish skin microbiomes (Boutin et al., 2014; Chiarello et al., 2015; Legrand et al., 2018; Krotman et al., 2020; Uren Webster et al., 2020), indicating that these phyla include symbiotic and commensal bacterial taxa that thrive on fish skin (Tarnecki et al., 2019). Another interesting finding was that three of the most abundant phyla that were present in both populations (Euryarchaeota, Deinococcus-Thermus, and Gemmatimonadetes; Figures 5A,B) were unique to fish skin microbiome samples (i.e., not detected in water), which could point to a specific association to fish microbiome. These three phyla have previously been reported in studies associated with fish microbiomes. For instance Gemmatimonadetes has been reported in the gut of Atlantic salmon (Salmo salar) (Dehler et al., 2017), and Euryarchaeota was detected in the skin microbiome of coral reef fishes (Chiarello et al., 2018). Deinococcus-Thermus have been detected in certain parts of the gut of the yellow grouper (Epinephelus awoara) (Egerton et al., 2018), but also on the skin of the common snook (Centropoumus undecimalis) (Tarnecki et al., 2021), and it is noteworthy that both of these species belong to the same order as perch, Perciformes.

The three ASVs present in at least 80% of the fish skin samples belonged to the families Burkholderiaceae and Rhizobiaceae, both of which are commonly associated with fish skin (Reinhart et al., 2019; Chiarello et al., 2020; Uren Webster et al., 2020). One ASV was annotated to species level, Variovorax paradoxus, a bacterium that can utilize many different organic compounds (Willems et al., 1991), and thus possibly a commensal microbiome member of perch that feed on the nutrient rich mucus (Shephard, 1994).

Taken together, these findings are in accordance with previous studies on fish skin microbiomes that also report on high variability among individuals in microbiome composition (Chiarello et al., 2015; Legrand et al., 2018, 2020; Uren Webster et al., 2020), and that fish skin microbiome composition is different from the surrounding bacterioplankton (Horsley, 1977; Wang et al., 2010; Stevens and Olson, 2013; Chiarello et al., 2018, 2019; Krotman et al., 2020; Sylvain et al., 2020; Uren Webster et al., 2020). One possible explanation for the different composition in fish skin microbiome compared to water, is that the fish skin offers a nutrient rich habitat for epibiotic bacteria in terms of the mucosal layer (Ángeles Esteban, 2012). The mucus constitutes a diversity of gel glycoproteins that have been suggested to provide many ecological niches for microbial organisms (Shephard, 1994; Wotton, 2004; Brown and Bythell, 2005; Chiarello et al., 2015, 2018; Sylvain et al., 2016).

Although we found statistical support for microbiome differences between populations, host individual accounted for a larger proportion of the total variation in microbiome composition than did host population. Individual-specific variation was evident also when considering the most abundant taxonomic groups; a few phyla were over- and underrepresented in the respective populations (e.g., Tenericutes and Patescibacteria, Figures 5A,B), however, this skewed abundance was an effect of host individual rather than population. A potential explanation for the relatively low differentiation between the two host populations in our study is that the host species is a strong determinant of the associated microbial community composition (Larsen et al., 2013; Stevens and Olson, 2013; Chiarello et al., 2018; Sylvain et al., 2020), at least relative to the effects of the environment (Östman et al., 2010). This was suggested by Stevens and Olson (2015) with the argument that the fish host-microbiome interaction creates niches that potentially make the skin surface less accessible for free-living microorganisms and possibly more resistant to fluctuations in abiotic factors (e.g., pH, temp, and nutrients) than bacterioplankton communities (Pinhassi et al., 2003; Martiny et al., 2006). It has previously been reported that population differences in skin microbiome compositions are correlated with genetic dissimilarities (Webster et al., 2018), and this could potentially apply to the pattern observed here since perch populations in the different study sites are genetically distinct (Olsson et al., 2011).

However, the high individual variation in microbiome composition, reported in our study, might be indicative of the fact that stochastic processes influence which microbes that colonize the skin of perch (Burns et al., 2016; Chiarello et al., 2019). Alternatively, the large variation of microbiomes among individual hosts may reflect individual differences in the genetic make-up of the immune system (Boutin et al., 2014; Malmstrøm et al., 2016; Webster et al., 2018), diet preferences (Chiarello et al., 2018; Uren Webster et al., 2020), and behaviors (Bolnick et al., 2003). Host individuals can be regarded as islands with different properties according to both intrinsic (e.g., host genetic and phenotypic variation) and extrinsic factors (e.g., external environmental conditions in the habitat). This means that their skin-associated microbiomes are exposed to constantly changing, and possibly contrasting, environmental conditions–both from the hosts and the environment that the host is exposed to. These combined features result in environmental heterogeneity which might have the potential to promote more species rich and diverse communities at higher level of biological organization (microbial communities), as have been shown for trophic levels in other systems (e.g., Koricheva and Hayes, 2018; Raffard et al., 2019).

In animal and plant ecosystems, heterogeneous or patchy environments are expected to promote species diversity (Johnson and Simberloff, 1974; Tews et al., 2004) and in this context, heterogeneity in microbiome composition within individuals can be of eco-evolutionary interest and importance (Anderson et al., 2006; Forsman and Wennersten, 2016; Yildirim et al., 2018). In line with this prediction, we hypothesized that different body parts constituted contrasting habitats and thus harbored different microbial communities (Costello et al., 2009; Chiarello et al., 2015). However, according to our results, the microbiomes on dorsal and ventral body sites within host individuals were not distinct from one another. This can either reflect similar microbiome assembly processes (Costello et al., 2009; Kraft et al., 2015), or that connectivity and dispersal of microbes between bodily regions within the host is high (MacArthur and Wilson, 1967; Miller et al., 2018).

Interestingly, we found that the within-individual variability of microbiome composition varied significantly among host individuals (see results from PERMDISP and Figure 2). This could reflect behavioral variation among host individuals. Vertical migration associated with foraging, thermoregulation and diel activity patterns expose the dorsal and ventral microbiomes to contrasting environmental conditions, and possibly also to different species pools of potential microbial colonizers (Bolnick et al., 2003; Nordahl et al., 2018, 2020). It can therefore be hypothesized that the difference between dorsal and ventral microbiomes should be more pronounced in individuals that engage in vertical migrations to a higher degree. Given that there is typically a large proportion of unexplained variation in microbiome composition across fish hosts (Falony et al., 2016; Miller et al., 2018), it is critical to evaluate the reliability and repeatability of microbiome samples taken from the same host individual. Knowledge about individual differences and measurement consistency and how they influence the partitioning of the total variance can inform sampling design, with potential to increase the reliability and to improve reproducibility of future studies (Voelkl et al., 2020). In our case, the high similarity of microbiomes sampled on the left and right side within perch individuals (Figure 6) suggests that the observed differences in variance between dorsal and ventral samples among individuals were not resulting from measurement error.

Previous attempts to sample, quantify and compare microbiomes between different functional parts of the fish host have not formally evaluated heterogeneity among body sites while accounting for individual identity using repeated samples from each body site (Chiarello et al., 2015; Lowrey et al., 2015; Legrand et al., 2018). The results from this study, based on repeated samples taken from the same individuals and body sites, thus provide novel insights on how variation among individuals in alpha diversity, and the degree of spatial heterogeneity within individuals, contribute to beta-diversity in fish skin microbiomes.