While the link between the gut microbiome and host behaviour is well established, how the microbiomes of other organs correlate with behaviour remains unclear. Additionally, behaviour–microbiome correlations are likely sex-specific because of sex differences in behaviour and physiology, but this is rarely tested. Here, we tested whether the skin microbiome of the Trinidadian guppy, Poecilia reticulata , predicts fish activity level and shoaling tendency in a sex-specific manner. High-throughput sequencing revealed that the bacterial community richness on the skin (Faith's phylogenetic diversity) was correlated with both behaviours differently between males and females. Females with richer skin-associated bacterial communities spent less time actively swimming. Activity level was significantly correlated with community membership (unweighted UniFrac), with the relative abundances of 16 bacterial taxa significantly negatively correlated with activity level. We found no association between skin microbiome and behaviours among male fish. This sex-specific relationship between the skin microbiome and host behaviour may indicate sex-specific physiological interactions with the skin microbiome. More broadly, sex specificity in host–microbiome interactions could give insight into the forces shaping the microbiome and its role in the evolutionary ecology of the host.

These methods are modified from the manuscript:

a)  Fish origin and maintenance

We used laboratory-bred descendants of guppies collected from the Caura River, Trinidad. Guppies were housed at densities of 1–2 fish per litre in 4.5-litre tanks in a recirculating system at 24 ± 1°C, on a 12 L:12 D lighting schedule and fed daily on flake and Artemia. All conditions were kept identical for the experiment duration, conducted in five experimental batches between August and November 2017.

b) Behavioural experiment and microbiome sampling

To assay behaviours, we used a central glass aquarium ('test aquarium'), flanked on both short sides by smaller aquaria ('stimulus aquaria'), all with 5cm water depth. The stimulus aquaria were lit from above by 32cm LED strip lights (350lm, 5W, 4000k, MeRox® Technics). The test aquarium had a black Perspex lid, to which was fastened a GoPro camera (Hero4 Black Edition; GoPro Inc. San Mateo, CA). We placed a test fish (of either sex: n=18 females; 19 males) in an acetate cylinder in the centre of the test tank for a 5-minute settling period. We then remotely set the GoPro to record, lifted the settling cylinder using a pulley, and recorded a 10-minute test period. For a subset of trials, we placed a shoal of three females in one of the stimulus tanks (side chosen haphazardly). The short sides of the test aquarium were one-way mirrors, such that test fish could see stimulus fish, but not vice versa. The experimental set-up was enclosed in a blackout curtain to ensure even lighting and prevent disturbance. The test aquarium was scrubbed with 70% ethanol, rinsed, and refilled between trials. We quantified "shoaling-tendency" as the proportion of trial time the test fish spent within the third of the test tank closest to the stimulus tank and "activity level" as the proportion of time the test fish actively swam throughout the 10 minutes. One of four observers, blind to the treatment, recorded the activity levels and shoaling levels from the videos using JWatcher™ 1.0 (www.jwatcher.ucla.edu).

c) Skin microbiome collection

The day after the behaviour trials, fish skin microbiome samples were collected. Fish were anaesthetised (0.02% tricaine methanesulfonate; PHARMAQ Ltd., Fordingbridge, UK), washed twice with two separate 20ml aliquots of sterile water, then swabbed the body surface with a sterile cotton swab for 30 seconds. Between fish, all glassware was sterilised with 70% ethanol. The swabs were stored at -20 ºC until DNA extraction.

d) Microbiome inventory

We used QIAamp PowerFecal DNA Isolation Kit (Qiagen, Hilden, Germany) to extract DNA, following the manufacturer's protocol. To control for contaminants, we included four blank samples in extractions and sequencing. Extracted DNA was shipped on dry ice to the DNA Services Facility at the University of Illinois at Chicago (Chicago, IL, USA). All library preparation, amplification, and sequencing were performed on the Illumina Miseq platform (V4 region of the bacterial 16S rRNA gene). We filtered the raw sequences for quality, removed chimeric reads, merged forward and reverse reads and assigned the remaining reads to ASVs (Amplicon Sequence Variant) using the DADA2 pipeline in QIIME2 2019.7. We trimmed each read to 220 base pairs, removed all singleton reads and used FastTree to build a phylogenetic tree, using the SILVA database classifier release 138 to assign taxonomy. We removed all non-bacterial reads (archaea, chloroplast, or mitochondria) and any ASVs detected in control blank samples prior to analysis.

d) Data analysis

First, we calculated Faith's phylogenetic distance (Faith's PD) as a metric for skin microbial alpha diversity for each sample, and compared Faith's PD between female and male guppies using a Kruskal–Wallis test using QIIME2.

To test for correlations between behaviour and microbiome alpha diversity, we used activity level or shoaling tendency as the response variable in two generalised linear mixed models (glmm) in the glmmTMB package, v. 1.1.2.3 (Beta error family and identity logit function within R v.4.1.0). As fixed main effects, we included test fish sex and length (controlling for sex differences), body condition (scaled mass index), Faith's PD, trial date and time of day, and whether there was a shoal present. We also included the two-way interaction between sex and Faith's PD. We used visreg and ggplot2 for plotting.

            To investigate correlates of skin microbiome bacterial community membership, we calculated the Unweighted UniFrac distance matrix in Qiime2, imported it into Primer v7, and built a distance-based linear model (DISTLM). This model included the same predictor variables as the glmm. We used permutation to assess statistical significance. We visualised the fitted DISTLM models in multidimensional space, using the distance-based redundancy analysis (dbRDA) in PRIMER v7 (see supplement).

Finally, to test for associations between individual bacterial taxa and guppy activity level and shoaling tendency, we used MaAsLin2. Within our models we added trial date as a fixed effect to control for any compositional differences caused by difference in experimental date. This programme first uses boosting in a univariate pre-screen to identify taxa (arcsine square-root transformed relative abundances) potentially associated with a focal variable and then uses linear mixed-effects models to test for statistically significant associations using the Benjamini–Hochberg false discovery rate method.

We have included .csv files containing the data, and a .pdf README file including the R code and output we used to analyse the data, plus more details on the variables in the datasheets.