Host-associated fungal communities are determined by host phylogeny and exhibit widespread associations with the bacterial microbiome

Interactions between hosts and their resident microbial communities are a fundamental component of fitness for both agents. Despite a recent proliferation of research on interactions between animals and their associated bacterial communities, comparative evidence from fungal communities is lacking, especially in natural populations. This disparity means knowledge of host-microbe interactions is biased towards the bacterial microbiome. Using samples from 49 species from eight metazoan classes, we demonstrate that the ecological distance between both fungal and bacterial components of the microbiome shift in tandem with host phylogenetic distance. Though so-called phylosymbiosis has been shown in bacterial communities, we extend previous knowledge by demonstrating that the magnitude of shifts in fungal and bacterial community structure across host phylogeny are correlated. These data are indicative of coordinated recruitment by hosts for specific suites of microbes, and potentially selection for bacterial-fungal interactions across a broad taxonomic range of host species. Using co-occurrence networks comprising both microbial groups, we illustrate that fungi form a critical component of microbial interaction networks, and that the strength and frequency of such interactions vary across host taxonomic groups. Collectively these data indicate fungal microbiomes may play a key role in host fitness and suggest an urgent need to study multiple agents of the animal microbiome to accurately determine the strength and ecological significance of host-microbe interactions.


INTRODUCTION
Multicellular organisms support diverse microbial communities that are critical for physiological functioning, immunity, development, evolution, behaviour, and even conservation [1][2][3] . Variability in demonstrated strong signals of phylosymbiosis across the broad host taxonomic range we test here.
Specifically, we predicted that patterns of phylosymbiosis within microbial kingdoms will also drive 76 significant positive covariance in patterns of microbial community structure between microbial 77 kingdoms within individual hosts, suggestive of evolutionary constraints that favour co-selection of 78 specific bacterial and fungal communities in tandem. Additionally, we used network analysis to identify 79 key bacteria-fungi interactions whilst quantifying variation in the frequency and strength of bacteria-80 fungi interaction networks across host taxonomic groups. Finally, we investigate the prediction that  were more likely to have higher bacterial relative to fungal diversity than Aves in our study organisms 98 (mean difference in probability 22.9% [1.6 -45.7%). Variation among species in this model explained models without sample preparation protocol effects, though this inflated the estimate of R 2 for all 132 taxonomic groupings (Table 2)

180
Our study represents the most wide-ranging evaluation of animal mycobiome composition, and its 181 covariation with the bacterial microbiome, undertaken to date. Our data provide novel evidence for illustrating frequent correlative links between fungal and bacterial taxa, whereby certain pairs of microbes from different kingdoms are much more likely to co-occur in the microbiome than expected 189 by chance. Taken together, these data provide novel evidence consistent with recruitment by animal 190 hosts for specific fungal and bacterial communities, which in turn may reflect selection for interactions 191 between bacteria and fungi critical for host physiology and health.

192
We found marked variation among host species in microbial community richness and 193 composition for both bacteria and fungi. Complementary analyses using mixed models and 194 phylogenetic models both detected a signal of host phylogeny in determining fungal and bacterial 195 microbiome diversity. Though our data suggest many species support a diverse assemblage of host-196 associated fungi, critically we show that i) bacterial diversity tends to be higher on average relative to 197 fungal diversity; and ii) there is no signal of positive covariance between fungal and bacterial richness 198 within species, suggesting more ASV-rich bacterial microbiomes are not consistently associated with 199 more ASV-rich mycobiomes. These patterns could arise because of competition for niche space

206
We detected strong phylosymbiosis for both fungi and bacteria across a broad host 207 phylogeny encompassing both vertebrate and invertebrate classes. This pattern was significantly 208 stronger in bacteria than for fungi. In both microbial kingdoms, the signal of phylosymbiosis 209 strengthened when aggregating microbial taxonomic assignments to family level, a phenomenon that 210 has previously been shown for bacterial communities 31 . That this pattern also occurs in fungi suggests 211 either that host recruitment is weaker at finer-scale taxonomies, or our ability to detect that signal is

220
In addition to microbe-specific patterns of phylosymbiosis, a key novel finding of our work is 221 that fungal and bacterial community composition correlate strongly across the host phylogeny. These 222 patterns are consistent with host recruitment for particular suites of fungal and bacterial taxa, which 223 may represent bacteria-fungi metabolic interactions beneficial to the host. Bacteria-fungi interactions 224 have previously been demonstrated for a handful of host species 8,9,17,33,34 , but here we show these are

241
The drivers of phylosymbiosis remain unclear, even for bacterial communities; is it a 242 phylogenetic signal indicative of host-microbiome coevolution, or simply a product of "ecological   (Table 1) and normalised to ~10 ng/ul. Samples were largely collated from previous studies 277 and/or those available from numerous researchers and as such, DNA extraction and storage 278 techniques were not standardised across species. We sequenced a median of 10 samples per species (range of 5 to 12; Table 1). MiSeq platform at the University of Salford. We ran the ITS rRNA library twice to increase sequencing 289 depth, and combined data within samples across these two runs in the data pre-processing stage.

290
We conducted all data processing and analysis in RStudio v1.2.1335 for R 52,53 (see 291 supplementary files for full code). We conducted amplicon sequence processing in DADA2 v1.5 54 for 292 both ITS rRNA and 16S rRNA amplicon data (see supplementary material for additional information).

293
After data processing, we obtained a median of 1425 reads per sample (range of 153 to 424,527) for 294 ITS rRNA libraries, and a median of 3273 reads per sample (range of 153 to 425,179) for 16S rRNA 295 libraries.

297
Host Phylogenetic Distances 298 As many of our host species lack genomic resources from which to construct a genome-based 299 phylogeny, we built a dated phylogeny of host species using TimeTree 55 . The phylogenetic tree 300 contained 42 species, of which 36 were directly represented in the TimeTree database. A further six 301 species had no direct match in TimeTree and so we used a congener as a substitute (Amietia, species based on shared branch length in the phylogeny using the 'cophenetic' function in the ape package 56 in R. We visualised and annotated the phylogeny using the R package ggtree 57 . fungal and bacterial kingdoms respectively. Alpha-diversity measures remained relatively stable within 316 a host species whether data were rarefied to 500, 1000, or 2500 reads, with similar patterns exhibited 317 between the two kingdoms (Figs. 1, S1, S2; see Supplementary Material for more details). To 318 visualise differences between microbial richness within species, we filtered the data to species with at 319 least two samples per microbial kingdom, giving a total of 41 species from six classes. For model 320 fitting, we filtered the data to only those samples with paired metrics of microbial richness for both 321 kingdoms (201 observations from 42 species). We fitted two models to these data. First, to quantify 322 relative differences in richness between bacteria and fungi within a sample, we used GLMMs in the 323 brms package, with i) Bernoulli errors and a logit link; ii) a binary response of '1' if bacterial richness 324 was higher than fungal richness, and '0' otherwise; and iii) 'Species' nested within 'Class' as random 325 intercepts. We did not include intermediate levels of taxonomy because replication at Order and 326 Family levels was low relative to Class. We did not use a phylogenetic mixed model as not all species 327 were represented in the TimeTree phylogeny. Second, to quantify absolute differences in microbial 328 richness, we fitted a bivariate response LMM with both fungal and bacterial richness values as a two-329 column response with Class as a fixed effect, and Species as a random intercept. For all models, we 330 used uninformative Cauchy priors for the random effects and Gaussian priors for fixed effects 331 coefficients. We assessed model adequacy using visual inspection of chains to assess mixing and stationarity properties, as well as posterior predictive checks using the 'pp_check' function in brms.
We calculated r 2 of models using the 'bayes_R2' function. We assessed the importance of terms 334 based on whether 95% credible intervals of the parameter estimates of interest crossed zero. We 335 used ggplot 60 , cowplot 61 and tidybayes 62 for raw data and model estimate plotting. To support these 336 analyses, we also used the R packages phylobase 63 and phylosignal 64 to estimate the phylogenetic 337 signal in patterns of alpha diversity for both bacteria and fungi, using both Inverse Simpson Index and 338 number of observed ASVs as outcome variables. We calculated Abouheif's Cmean for each diversity-339 microbe combination and corrected p values for multiple testing using Benjamini-Hochberg correction.

340
To identify taxonomic differences in microbiome and mycobiome composition between host 341 species, we used centred-log-ratio (CLR) transformation in the microbiome 65 package to normalise 342 microbial abundance data, which obviates the need to lose data through rarefying 66 . To visualise

356
To test the hypothesis that inter-individual differences in microbial community composition 357 were preserved between microbial kingdoms, we performed Procrustes rotation of the two PCA 358 ordinations for bacterial and fungal abundance matrices, respectively (n = 277 paired samples from 359 46 species). We also repeated this analysis with ASVs agglomerated into progressively higher each kingdom at each iteration, we randomly sampled 90% of the data and recalculated the correlation metric. We repeated this process 999 times to build a distribution of correlation values at 364 each taxonomic grouping.

365
To examine the hypothesis that inter-individual distance in microbial community composition 366 varies in concert with inter-host phylogenetic distance, we performed a Procrustes rotation on the 367 paired matrix of microbial distance (Euclidean distance of CLR-transformed abundances) and patristic 368 distance from the phylogenetic tree. We also repeated the analysis but binning the microbial data to 369 family level. As above, we conducted a bootstrap resampling procedure, selecting 90% of the data at 370 random and recalculating the correlation, for a total of 1000 permutations. This allowed us to test for 371 significant differences in strength of correlation within kingdoms across taxonomic grouping levels,

372
and across kingdoms within a particular taxonomic grouping.

373
To determine the effect of diet on bacterial and fungal community composition, we used only 374 samples from the bird and mammal species and agglomerated the data for each host species using 375 the merge_samples function in phyloseq 68 . This gave us an representative microbiome for each host 376 species, which we rarefied to the lowest number of reads for each combination of kingdom and host 377 taxon (2,916 -9,160 reads; bacterial read counts were low for lesser horseshoe bats and so this 378 species was removed from this analysis) and extracted Euclidean distance matrices for each. We

383
To identify potential relationships between fungal and bacterial communities, we conducted 384 two analyses; 1) We used the R package SpiecEasi 70 to identify correlations between unrarefied, 385 CLR-transformed ASVs abundances at the host class level (with Invertebrates grouped), and 2) we 386 used co-occurrence analysis at the species level, by rarefying the bacterial and fungal data sets to 387 500 reads each, and agglomerated these to family level, resulting in 117 bacterial groups and 110 388 fungal groups. We then merged the phyloseq objects for bacterial and fungal communities for each constructing a Spearman's correlation coefficient matrix in the bioDist package 71,72 . We visualised the class-level microbial networks comprising both positive and negative interactions using the 395 modularity() function after greedy clustering implemented in the igraph package. We resampled 90% 396 of nodes in each network 1000 times to build distributions of modularity with which to quantify 397 differences among animal classes. We used binomial GLM to test the hypothesis that the proportion