Variation in the microbiome in response to parasite, bat host, and environment
Ordination of microbiomes and PERMANOVA indicated that the parasite (i.e., parasite family and species), the host bat (i.e., bat family, bat sex, and bat individual), and the environment (i.e., region and sampling site) significantly contributed to bat fly microbiome variation (Table 2; Figure 2B and Figure S6). Other variables significantly contributed to microbiome community differentiation (i.e., bat feeding guild, bat species, protection status of sampling site, habitat patch area, and isolation), but violated PERMANOVA’s assumption of homoscedasticity. Parasite species, parasite family, bat feeding guild, bat family, bat species, and sampling site had the largest effect sizes, however many of these variables are correlated with each other.
Sequential sum-of-squares with free permutation indicated that parasite species significantly impacted microbiome community structure, consistent with single variable PERMANOVA (Table 3). Habitat patch protection status and region also significantly contributed to microbiome variation. In the four most well-sampled species, none of the test variables significantly explained microbiome variation without violating PERMANOVA assumptions (Table S2; Figure S7).
Bacterial taxon richness across habitat patches
While the three sampled sites within REGUA had the highest ASV richness, there was no pattern of decreasing ASV richness with decreasing area, increasing isolation, or increasing distance from a source (Spearman; Area: rho=0.0613, p-value=0.4031; Isolation: rho=0.0020, p-value=0.9787; Distance to Source: rho=-0.0564, p-value=0.4421; Figure 3). Median bacterial ASV richness fell between 6 and 11 for each sampled parasite individual, but the range of ASV richness per fragment varied dramatically. There was no significant difference between sampling sites in mean ASV richness (Kruskal-Wallis chi-squared statistic=17.856, p-value=0.1201).
Bacterial association networks
The connectance and betweenness centrality of habitat patch networks varied with habitat patch effects, but modularity did not. Network connectance was significantly lower in REGUA patches than in networks constructed for sites outside of REGUA (Figure 5B; Wilcoxon Rank Sum test, p-value=0.012). Betweenness centrality was higher in networks for patch F4, F2, and F1 (i.e., smallest patches), corresponding with their band-like network structure (Figure 4), and centrality was significantly different between networks from large and small patches (Figure 5c; Kruskal-Wallis chi-squared statistic = 300.23, p-value < 2.2 x 10-16). Raw and Z-score modularity of habitat patch networks did not significantly differ between REGUA and non-REGUA sites, but tended to be lower in patches outside REGUA (Figure S8). Z-score modularity using the measured networks for standardization did not control for network size and shape, and mimicked the pattern exhibited by raw modularity. Modularity measures standardized by null distributions did not support patterns of decreasing modularity with decreasing patch area (Figure S8).
In habitat patch networks reconstructed for the three REGUA sites and large patches, low-abundance bacteria had higher betweenness centrality than other bacteria in these networks, acting to connect graph neighborhoods (Figure 4). As habitat patch networks became less modular with decreasing habitat patch area, nodes with high betweenness centrality were typically low-abundance bacteria and bacteria in the genera Arsenophonus (including Candidatus Phlomobacter), Wolbachia, and Bartonella.
While there was no impact of sample size (i.e., number of parasite individuals) on ASV richness in each network, a greater number of samples may allow detection of more edges between nodes and change the size of a network. Sample size was lowest in small fragments and highest in large fragments, with the exception of F4 which has intermediate area, isolation, and distance to source measurements, but supports a high diversity of parasites. We did not detect a significant correlation between sample size and network modularity (Spearman rho=0.3945, p-value=0.229). Using the subsampling scheme within SpiecEasi, we did not detect a correlation of sample size and proportion of high-confidence network edges (Spearman rho=0.0275, p-value=0.936).
Using ordinations of graphlet orbits as a size-independent comparison of networks [84, 88, 89], habitat patch networks vary with habitat area (Figure 5D) and species-specific networks are distinct within and outside of REGUA (Figure S8). Principal coordinate axis 2 corresponded well to decreasing habitat area, with large fragments positioned higher on the access and gradually decreasing in patch area lower on the axis. Principal coordinate axis 1 primarily illustrated variation in networks from patches F8 and F5. This variation does not correspond with environment, parasite, or bat variables. K-means clustering with 3 groups separated large habitat patch networks (REGUA networks, F10, F9, F7) from small patch networks (F4, F2, F1), and F8 and F5 formed a unique cluster. Species-specific networks also indicated distinctions between within-REGUA and outside-REGUA networks. All within-REGUA networks occupied unique ordination space from outside-REGUA networks within species.
Discussion
Our research builds upon previous evidence that the environment influences microbiome composition in addition to host factors [90–93]. However, previous studies have been primarily conducted on the microbiomes of free-living hosts or microbiomes within a single host species, preventing the examination of how the environment might impact communities of interacting macro- and microorganisms. Contrary to expectations, habitat loss did not lead to a decrease in bacterial diversity. Instead, habitat loss impacted bacterial community structure, which includes diversity and relative abundance. The interactions of bat flies with the broader environment are filtered through their obligate associations with host bats, yet the signal of environmental change is also detected in the composition of bat fly microbiomes. This indicates that environmental degradation may have cascading consequences through hierarchical communities of interacting organisms.
Parasite and environment as drivers of microbiome composition
Parasite species identity is the strongest predictor of microbiome composition (Tables 2 and 3; Figure 2B). Specificity of arthropod microbiomes has been previously found in tsetse flies (Glossinidae; [94], which are also members of the Hippoboscoidea. That microbiome composition shows such a strong signal of parasite species identity indicates that either more of the microbiome is maternally inherited than previously understood, or that even non-maternally inherited bacteria may be maintained through life history traits (e.g., host bat associations, microclimatic preference; horizontal transmission of bacteria).
Habitat patch area and protection status, but not degree of isolation or distance to source, have a measurable effect on the microbiome of bat flies, but less explanatory power than parasite species (Tables 2 and 3; Figures 3, 4, 5). While bacterial ASV richness does not vary following expectations of island biogeography theory (Figure 3), examining both relative abundance and diversity of bacteria using PERMANOVA and bacterial association networks provided a clear statistical signal that habitat patch area (measured as area and protection status) is correlated with microbiome composition (Table 3; Figure 5). Other measures like isolation and distance to source had no impact on bat, parasite, or microbiome communities, likely because bats can easily move between these patches. However, even when habitat patches are proximal, fragmented landscapes in the Atlantic Forest have been previously shown to significantly impact bat-plant and bat-ectoparasite interaction networks [95].
Even though changes in the bat and parasite communities explain most of the variation in microbiome association networks, the environment may also directly impact bat fly microbiome communities. The structure of habitat patch networks covaried with habitat area and parasite species-specific networks were consistently different within REGUA compared to outside REGUA. In addition, parasite and bat species diversity did not always correspond to a specific network structure. REGUA and fragments F10 (the largest fragment outside of REGUA, 9 bat fly species) and F4 (8 bat fly species) have the greatest parasite richness, but F4 is intermediate in area and isolation. Networks from patch F4 consistently cluster with other networks from small patches (Figure 5C,D), despite having similar parasite species richness to large fragments. That patch F4 has high parasite species richness but similar bacterial association networks to smaller fragments (i.e., lower parasite species richness), suggests that bat and parasite community composition does not solely explain microbiome variation. The environment may also directly drive composition of parasite-associated microbiomes.
Implications of changes in bacterial network structure in response to environment
Bacteria with high betweenness centrality may act as hub species that maintain the stability of a network [19, 92, 93]. In small fragments, Arsenophonus, Wolbachia, and Bartonella had high betweenness centrality, but these bacterial taxa were less central to the networks from large fragments despite maintaining high relative abundance in flies at these sites. Decreasing betweenness centrality may be indicative of changing interactions between bacteria in response to environmental perturbations. Bacteria of blood-feeding insects play an important role in vector competence in insects [19, 96, 97]. For example, in tsetse flies, primary bacterial endosymbionts in the genus Wigglesworthia impede the invasion of trypanosome parasites by assisting host defenses and subsequently decrease the competence of tsetse flies to vector these harmful parasites to downstream hosts including humans [98]. As bat flies are important arthropod vectors of bat pathogens [99], changes in the structure of their microbiomes in response to habitat loss may have implications for the disease ecology of arthropod vectored pathogens in bats.
Modules may delimit groups of bacteria with specific functional specializations and/or groups that respond in similar ways to environmental variables [99, 100]. Higher modularity may protect a community of free-living organisms from invading pathogens, because a pathogen would be isolated to one module within the community (i.e., diversity-stability debate) [100, 101]. This hypothesis may be applicable to bacterial networks if pathogens are limited in transmission by direct competition with endogenous bacteria. However, high modularity in bacterial association networks may also reflect the absence of microbiome-mediated host defenses against pathogen invasion. If we revisit the tsetse fly example, Wigglesworthia and trypanosome parasites would share an edge in a microbe association network because Wigglesworthia abundance and prevalence is associated with low trypanosome abundance and prevalence mediated by host defenses. The absence of this interaction in a network might indicate that Wigglesworthia-induced host response to trypanosome invasion is impaired by other microbes in the network (e.g., priority effects) or other aspects of host health. Trypanosome transmission would be less regulated in this instance, but modularity would be high as long as trypanosome abundance and prevalence did not lead to changes in the abundance or prevalence of other microbes. In the case of bat flies, the consequences of more isolated modules for the functionality and stability of the microbiome are unclear and merit future investigation.
The missing primary symbionts of Brazilian nycteribiid flies
Arsenophonus has been previously identified as a primary symbiont of streblid bat flies and the closely-related “Candidatus Aschnera chinzeii” is the hypothesized primary symbiont of some nycteribiid bat flies [49, 50, 102–104]. The high relative abundance and prevalence of Arsenophonus in all streblid bat flies sampled in our study is consistent with the hypothesis that Arsenophonus acts as the primary symbiont in the family Streblidae. However, it is unlikely that Arsenophonus or “Candidatus Aschnera chinzeii” are acting as the primary symbionts of the nycteribiid species that we sampled (Figure 2A). Previous studies that have identified “Candidatus Aschnera chinzeii” as the primary symbiont in nycteribiid flies examined species that are geographically limited to Africa, Asia, Europe, and Oceania [102, 103]. Species within Basilia, the only globally distributed bat fly genus [105, 106], have varying symbiont associations [49, 50, 103, 104]. Only one previous study has examined symbionts of Basilia from The Americas [107]. This study detected an Arsenophonus variant in two individuals of Basilia boardmani (Nycteribiinae; restricted to North America) that was distinct from “Candidatus Aschnera chinzeii” and from Arsenophonus detected in other nyteribiid species. In the Basilia species sampled from The Atlantic Forest, of which 3 are limited to South America (B. andersoni, B. juquiensis, B. lindolphoi) and 1 is found in North and South America (B. ferruginea), we detect no “Candidatus Aschnera chinzeii” and low relative abundance or an absence of Arsenophonus (including “Candidatus Phlomobacter”). These findings in combination with previous studies suggest that neither “Candidatus Aschnera chinzeii” nor Arsenophonus act as the primary symbiont of the sampled Basilia species. It may be that Wolbachia and Bartonella both have high relative abundance in nycteribiid bat flies because they are acting as primary or facultative symbionts, or because they are acting as a reproductive parasite and pathogen, respectively [107–113]. If the latter is true, the primary symbiont may be one of the less abundant bacteria or another microbe not detected by 16S rRNA sequencing (e.g. fungi; [114].
Summation
Understanding of community structure and function must be extended beyond an exclusively macroorganismal view to include the layers of microorganisms that also participate in defining an ecological community. By examining the hierarchical interactions between bats, bat flies, and bat fly bacterial microbiomes within a largely deforested landscape, we attempted to more accurately characterize the consequences of environmental change for a wildlife community. Reduced fragment area led to generally less diverse and less abundant bat fly communities, which led to less modular bacterial association networks, but without a decrease in bacterial ASV richness. Our data highlight the importance of considering ecological responses of microorganismal taxa. Future research is needed in order to capture the full impact of continuing deforestation and habitat change on ecological communities.
Acknowledgements
Thank you to REGUA staff for their assistance in the field and Dr. Marcelo Weksler and Museu Nacional/UFRJ staff for logistical help. This research was funded by Richard Gilder Graduate School Dissertation Research Fellowships to KAS and MRI, and a CAPES “science without borders” Fellowship to TSMT with additional funding from Queen Mary University of London to ELC and from The American Museum of Natural History to SLP. This research was also enabled by funding to ELC from the Natural Sciences and Engineering Research Council of Canada through the Discovery Grants Program, support from the Government of Canada's New Frontiers in Research Fund (NFRFT-2020-0073) and by support from Genome Canada and Ontario Genomics to BIOSCAN–Canada through the Large Scale Applied Research Program. All field work was conducted under Instituto Brasileiro do Meio Ambiente e dos Recursos Naturais Renováveis permit 19037-1.
Conflict of Interest Statement
The authors have no competing interests to declare.
Data Accessibility Statement
Raw 16S rRNA reads are available under NCBI SRA accessions SRX13352735-13352962 and BioProject PRJNA786937, COI sequences of bat flies are available on GenBank (OL847352-OL847639), metadata and all code used for analyses is available on Github (DOI: 10.5281/zenodo.6796123).
Literature Cited
1. Alroy J. Effects of habitat disturbance on tropical forest biodiversity. Proc Natl Acad Sci U S A 2017; 114: 6056–6061.
2. Gatti LV, Basso LS, Miller JB, Gloor M, Gatti Domingues L, Cassol HLG, et al. Amazonia as a carbon source linked to deforestation and climate change. Nature 2021; 595: 388–393.
3. Ellwanger JH, Kulmann-Leal B, Kaminski VL, Valverde-Villegas JM, Veiga ABGDA, Spilki FR, et al. Beyond diversity loss and climate change: Impacts of Amazon deforestation on infectious diseases and public health. An Acad Bras Cienc 2020; 92: e20191375.
4. Morand S, Lajaunie C. Outbreaks of vector-borne and zoonotic diseases are associated with changes in forest cover and oil palm expansion at global scale. Front Vet Sci 2021; 8: 661063.
5. Keenan RJ, Reams GA, Achard F, de Freitas JV, Grainger A, Lindquist E. Dynamics of global forest area: Results from the FAO Global Forest Resources Assessment 2015. For Ecol Manage 2015; 352: 9–20.
6. Rezende CL, Scarano FR, Assad ED, Joly CA, Metzger JP, Strassburg BBN, et al. From hotspot to hopespot: An opportunity for the Brazilian Atlantic Forest. Perspectives in Ecology and Conservation 2018; 16: 208–214.
7. Yarwood SA. The role of wetland microorganisms in plant-litter decomposition and soil organic matter formation: a critical review. FEMS Microbiol Ecol 2018; 94.
8. Kock RA, Orynbayev M, Robinson S, Zuther S, Singh NJ, Beauvais W, et al. Saigas on the brink: Multidisciplinary analysis of the factors influencing mass mortality events. Sci Adv 2018; 4: eaao2314.
9. Murdock CC, Blanford S, Hughes GL, Rasgon JL, Thomas MB. Temperature alters Plasmodium blocking by Wolbachia. Sci Rep 2014; 4: 3932.
10. MacArthur RH, Wilson EO. An equilibrium theory of insular zoogeography. Evolution 1963; 17: 373–387.
11. Krasnov BR, Shenbrot GI, Medvedev SG. Host–habitat relations as an important determinant of spatial distribution of flea assemblages (Siphonaptera) on rodents in the Negev Desert. Parasitology 1997.
12. Poulin R. Are there general laws in parasite ecology? Parasitology 2007.
13. Speer KA, Dheilly NM, Perkins SL. Microbiomes are integral to conservation of parasitic arthropods. Biol Conserv 2020; 108695.
14. Bell T, Ager D, Song J-I, Newman JA, Thompson IP, Lilley AK, et al. Larger islands house more bacterial taxa. Science 2005; 308: 1884.
15. Zinger L, Boetius A, Ramette A. Bacterial taxa-area and distance-decay relationships in marine environments. Mol Ecol 2014; 23: 954–964.
16. Martiny JBH, Bohannan BJM, Brown JH, Colwell RK, Fuhrman JA, Green JL, et al. Microbial biogeography: putting microorganisms on the map. Nat Rev Microbiol 2006; 4: 102–112.
17. Carbonero F, Oakley BB, Purdy KJ. Metabolic flexibility as a major predictor of spatial distribution in microbial communities. PLoS One 2014; 9: e85105.
18. van der Gast CJ. Microbial biogeography: the end of the ubiquitous dispersal hypothesis? Environ Microbiol . 2015. , 17: 544–546
19. Weiss B, Aksoy S. Microbiome influences on insect host vector competence. Trends Parasitol 2011; 27: 514–522.
20. Gupta A, Nair S. Dynamics of insect-microbiome interaction influence host and microbial symbiont. Front Microbiol 2020; 11: 1357.
21. Dick CW, Dittmar K. Parasitic bat Flies (Diptera: Streblidae and Nycteribiidae): Host specificity and potential as vectors. In: Klimpel S, Mehlhorn H (eds). Bats (Chiroptera) as Vectors of Diseases and Parasites. 2014. Springer, Berlin, Heidelberg, pp 131–155.
22. Speer KA, Luetke E, Bush E, Sheth B, Gerace A, Quicksall Z, et al. A fly on the cave wall: Parasite genetics reveal fine-scale dispersal patterns of bats. 2019.
23. Patterson BD, Dick CW, Dittmar K. Sex biases in parasitism of neotropical bats by bat flies (Diptera: Streblidae). J Trop Ecol 2008; 24: 387–396.
24. Hiller T, Brändel SD, Honner B, Page RA, Tschapka M. Parasitization of bats by bat flies (Streblidae) in fragmented habitats. Biotropica 2020; 72: 617.
25. Kikuchi Y, Tada A, Musolin DL, Hari N, Hosokawa T, Fujisaki K, et al. Collapse of insect gut symbiosis under simulated climate change. MBio 2016; 7.
26. Thapa S, Zhang Y, Allen MS. Effects of temperature on bacterial microbiome composition in Ixodes scapularis ticks. Microbiologyopen 2019; 8: e00719.
27. Souto Martins Teixeira T. Bats in a fragmented world. 2019. Queen Mary University of London.
28. Emmons L, Feer F. Neotropical rainforest mammals: a field guide. 1997. sidalc.net.
29. Reis NR, Fregonezi MN, Peracchi AL, Shibatta OA. Morcegos do Brasil: guia de campo. 2013. Technical Books Editora.
30. Sikes RS, Care A, of Mammalogists UC of TAS. 2016 Guidelines of the American Society of Mammalogists for the use of wild mammals in research and education. J Mammal 2016; 97: 663–688.
31. Wenzel RL. The streblid batflies of Venezuela (Diptera: Streblidae). Brigham Young University Science Bulletin, Biological Series 1976; 20: 1.
32. Graciolli G, de Carvalho CJB. Moscas ectoparasitas (Diptera, Hippoboscoidea) de morcegos (Mammalia, Chiroptera) do Estado do Paraná. 11. Streblidae. Chave pictórica para gêneros e espécies 1. RevIa bras Zool 2001; 18: 907–960.
33. Graciolli G, de Carvalho CJB. Moscas ectoparasitas (Diptera, Hippoboscoidea, Nycteribiidae) de morcegos (Mammalia, Chiroptera) do Estado do Paraná, Brasil. I. Basilia, taxonomia e chave pictórica para as espécies 1. RevIa bras Zool 2001; 18: 33–49.
34. Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol 1994; 3: 294–299.
35. Hebert PDN, Cywinska A, Ball SL, deWaard JR. Biological identifications through DNA barcodes. Proceedings of the Royal Society of London Series B: Biological Sciences 2003; 270: 313–321.
36. Gustafson EJ, Parker GR. Relationships between landcover proportion and indices of landscape spatial pattern. Landsc Ecol 1992; 7: 101–110.
37. McGarigal K, Cushman SA, Neel MC, Ene E. FRAGSTATS: spatial pattern analysis program for categorical maps. 2002.
38. Gilbert JA, Meyer F, Antonopoulos D, Balaji P, Brown CT, Brown CT, et al. Meeting report: the terabase metagenomics workshop and the vision of an Earth microbiome project. Stand Genomic Sci 2010; 3: 243–248.
39. Gilbert JA, Jansson JK, Knight R. The Earth Microbiome project: successes and aspirations. BMC Biol 2014; 12: 69.
40. Apprill A, McNally S, Parsons R, Weber L. Minor revision to V4 region SSU rRNA 806R gene primer greatly increases detection of SAR11 bacterioplankton. Aquat Microb Ecol 2015; 75: 129–137.
41. Parada AE, Needham DM, Fuhrman JA. Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ Microbiol 2016; 18: 1403–1414.
42. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 2011; 17: 10–12.
43. Katoh K, Misawa K, Kuma K-I, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res 2002; 30: 3059–3066.
44. Katoh K, Kuma K-I, Toh H, Miyata T. MAFFT version 5: improvement in accuracy of multiple sequence alignment. Nucleic Acids Res 2005; 33: 511–518.
45. Katoh K, Toh H. PartTree: an algorithm to build an approximate tree from a large number of unaligned sequences. Bioinformatics 2007; 23: 372–374.
46. Price MN, Dehal PS, Arkin AP. FastTree 2--approximately maximum-likelihood trees for large alignments. PLoS One 2010; 5: e9490.
47. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods 2016; 13: 581–583.
48. DeSantis TZ, Hugenholtz P, Larsen N, Rojas M, Brodie EL, Keller K, et al. Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl Environ Microbiol 2006; 72: 5069–5072.
49. Hosokawa T, Nikoh N, Koga R, Satô M, Tanahashi M, Meng X-Y, et al. Reductive genome evolution, host–symbiont co-speciation and uterine transmission of endosymbiotic bacteria in bat flies. ISME J 2012; 6: 577–587.
50. Duron O, Schneppat UE, Berthomieu A, Goodman SM, Droz B, Paupy C, et al. Origin, acquisition and diversification of heritable bacterial endosymbionts in louse flies and bat flies. Mol Ecol 2014; 23: 2105–2117.
51. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 2013; 41: D590-6.
52. Yilmaz P, Parfrey LW, Yarza P, Gerken J, Pruesse E, Quast C, et al. The SILVA and “All-species Living Tree Project (LTP)” taxonomic frameworks. Nucleic Acids Res 2014; 42: D643-8.
53. Glöckner FO, Yilmaz P, Quast C, Gerken J, Beccati A, Ciuprina A, et al. 25 years of serving the community with ribosomal RNA gene reference databases and tools. J Biotechnol 2017; 261: 169–176.
54. Nováková E, Hypsa V, Moran NA. Arsenophonus, an emerging clade of intracellular symbionts with a broad host distribution. BMC Microbiol 2009; 9: 143.
55. Bressan A, Terlizzi F, Credi R. Independent origins of vectored plant pathogenic bacteria from arthropod-associated Arsenophonus endosymbionts. Microb Ecol 2012; 63: 628–638.
56. Salter SJ, Cox MJ, Turek EM, Calus ST, Cookson WO, Moffatt MF, et al. Reagent and laboratory contamination can critically impact sequence-based microbiome analyses. BMC Biol 2014; 12: 87.
57. Weiss S, Amir A, Hyde ER, Metcalf JL, Song SJ, Knight R. Tracking down the sources of experimental contamination in microbiome studies. Genome Biol 2014; 15: 564.
58. Eisenhofer R, Minich JJ, Marotz C, Cooper A, Knight R, Weyrich LS. Contamination in low microbial biomass microbiome studies: Issues and recommendations. Trends Microbiol 2019; 27: 105–117.
59. Davis NM, Proctor DM, Holmes SP, Relman DA, Callahan BJ. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome 2018; 6: 226.
60. Bokulich NA, Subramanian S, Faith JJ, Gevers D, Gordon JI, Knight R, et al. Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing. Nat Methods 2013; 10: 57–59.
61. Alberdi A, Aizpurua O, Gilbert MTP, Bohmann K. Scrutinizing key steps for reliable metabarcoding of environmental samples. Methods Ecol Evol 2018; 9: 134–147.
62. McMurdie PJ, Holmes S. Phyloseq: a bioconductor package for handling and analysis of high-throughput phylogenetic sequence data. Pac Symp Biocomput 2012; 235–246.
63. McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 2013; 8: e61217.
64. Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: community ecology package. R package version 2.4--4. 2017. 2018.
65. Wickham H. ggplot2: Elegant Graphics for Data Analysis. 2016. Springer.
66. Fernandes AD, Reid JN, Macklaim JM, McMurrough TA, Edgell DR, Gloor GB. Unifying the analysis of high-throughput sequencing datasets: characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis. Microbiome 2014; 2: 15.
67. Gloor GB, Reid G. Compositional analysis: a valid approach to analyze microbiome high-throughput sequencing data. Can J Microbiol 2016; 62: 692–703.
68. Tsilimigras MCB, Fodor AA. Compositional data analysis of the microbiome: fundamentals, tools, and challenges. Ann Epidemiol 2016; 26: 330–335.
69. Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. Microbiome datasets are compositional: And this is not optional. Front Microbiol 2017; 8: 2224.
70. Silverman JD, Washburne AD, Mukherjee S, David LA. A phylogenetic transform enhances analysis of compositional microbiota data. Elife 2017; 6.
71. Xia Y, Sun J. Hypothesis testing and statistical analysis of microbiome. Genes Dis 2017; 4: 138–148.
72. Anderson MJ, Walsh DCI. PERMANOVA, ANOSIM, and the Mantel test in the face of heterogeneous dispersions: What null hypothesis are you testing? Ecol Monogr 2013; 83: 557–574.
73. Anderson MJ. Permutational Multivariate Analysis of Variance (PERMANOVA). In: Balakrishnan N, Colton T, Everitt B, Piegorsch W, Ruggeri F, Teugels JL (eds). Wiley StatsRef: Statistics Reference Online. 2014. John Wiley & Sons, Ltd, Chichester, UK, pp 1–15.
74. Kurtz ZD, Müller CL, Miraldi ER, Littman DR, Blaser MJ, Bonneau RA. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput Biol 2015; 11: e1004226.
75. Liu H, Roeder K, Wasserman L. Stability Approach to Regularization Selection (StARS) for High Dimensional Graphical Models. Adv Neural Inf Process Syst 2010; 24: 1432–1440.
76. Röttjers L, Faust K. From hairballs to hypotheses-biological insights from microbial networks. FEMS Microbiol Rev 2018; 42: 761–780.
77. Freeman LC. Centrality in social networks conceptual clarification. Soc Networks 1978; 1: 215–239.
78. Brandes U. A faster algorithm for betweenness centrality. J Math Sociol 2001; 25: 163–177.
79. Csardi G, Nepusz T, Others. The igraph software package for complex network research. InterJournal, Complex Systems 2006; 1695: 1–9.
80. Newman MEJ. Finding community structure in networks using the eigenvectors of matrices. Phys Rev E Stat Nonlin Soft Matter Phys 2006; 74: 036104.
81. Delmas E, Besson M, Brice M-H, Burkle LA, Dalla Riva GV, Fortin M-J, et al. Analysing ecological networks of species interactions: Analyzing ecological networks. Biol Rev 2019; 94: 16–36.
82. Fortunato S, Hric D. Community detection in networks: A user guide. arXiv [physics.soc-ph] . 2016.
83. Singh A, Humphries MD. Finding communities in sparse networks. Sci Rep 2015; 5: 8828.
84. Yaveroğlu ÖN, Malod-Dognin N, Davis D, Levnajic Z, Janjic V, Karapandza R, et al. Revealing the hidden language of complex networks. Sci Rep 2014; 4: 4547.
85. Przulj N. Biological network comparison using graphlet degree distribution. Bioinformatics 2007; 23: e177-83.
86. Hočevar T, Demšar J, Others. Computation of graphlet orbits for nodes and edges in sparse graphs. J Stat Softw 2016; 71.
87. Müller CL, Bonneau R, Kurtz Z. Generalized stability approach for regularized graphical models. arXiv [statME] . 2016.
88. Mahana D, Trent CM, Kurtz ZD, Bokulich NA, Battaglia T, Chung J, et al. Antibiotic perturbation of the murine gut microbiome enhances the adiposity, insulin resistance, and liver disease associated with high-fat diet. Genome Med 2016; 8: 48.
89. Ruiz VE, Battaglia T, Kurtz ZD, Bijnens L, Ou A, Engstrand I, et al. A single early-in-life macrolide course has lasting effects on murine microbial network topology and immunity. Nat Commun 2017; 8: 518.
90. Amato KR, Yeoman CJ, Kent A, Righini N, Carbonero F, Estrada A, et al. Habitat degradation impacts black howler monkey (Alouatta pigra) gastrointestinal microbiomes. ISME J 2013; 7: 1344–1353.
91. Avena CV, Parfrey LW, Leff JW, Archer HM, Frick WF, Langwig KE, et al. Deconstructing the bat skin microbiome: Influences of the host and the environment. Front Microbiol 2016; 7: 1753.
92. Becker CG, Longo AV, Haddad CFB, Zamudio KR. Land cover and forest connectivity alter the interactions among host, pathogen and skin microbiome. Proc Biol Sci 2017; 284.
93. Ingala MR, Becker DJ, Bak Holm J, Kristiansen K, Simmons NB. Habitat fragmentation is associated with dietary shifts and microbiota variability in common vampire bats. Ecol Evol 2019; 9: 6508–6523.
94. Aksoy E, Telleria EL, Echodu R, Wu Y, Okedi LM, Weiss BL, et al. Analysis of multiple tsetse fly populations in Uganda reveals limited diversity and species-specific gut microbiota. Appl Environ Microbiol 2014; 80: 4301–4312.
95. Mello RM, Laurindo RS, Silva LC, Pyles MV, Mancini MCS, Dáttilo W, et al. Landscape configuration and composition shape mutualistic and antagonistic interactions among plants, bats, and ectoparasites in human-dominated tropical rainforests. Acta Oecol 2021; 112: 103769.
96. Cirimotich CM, Ramirez JL, Dimopoulos G. Native microbiota shape insect vector competence for human pathogens. Cell Host Microbe 2011; 10: 307–310.
97. Sassera D, Epis S, Pajoro M, Bandi C. Microbial symbiosis and the control of vector-borne pathogens in tsetse flies, human lice, and triatomine bugs. Pathog Glob Health 2013; 107: 285–292.
98. Weiss BL, Wang J, Maltz MA, Wu Y, Aksoy S. Trypanosome infection establishment in the tsetse fly gut is influenced by microbiome-regulated host immune barriers. PLoS Pathog 2013; 9: e1003318.
99. Obame-Nkoghe J, Rahola N, Bourgarel M, Yangari P, Prugnolle F, Maganga GD, et al. Bat flies (Diptera: Nycteribiidae and Streblidae) infesting cave-dwelling bats in Gabon: diversity, dynamics and potential role in Polychromophilus melanipherus transmission. Parasit Vectors 2016; 9: 333.
100. Krause AE, Frank KA, Mason DM, Ulanowicz RE, Taylor WW. Compartments revealed in food-web structure. Nature 2003; 426: 282–285.
101. Stouffer DB, Bascompte J. Understanding food-web persistence from local to global scales. Ecol Lett 2010; 13: 154–161.
102. Trowbridge RE, Dittmar K, Whiting MF. Identification and phylogenetic analysis of Arsenophonus- and Photorhabdus-type bacteria from adult Hippoboscidae and Streblidae (Hippoboscoidea). J Invertebr Pathol 2006; 91: 64–68.
103. Morse SF, Bush SE, Patterson BD, Dick CW, Gruwell ME, Dittmar K. Evolution, multiple acquisition, and localization of endosymbionts in bat flies (Diptera: Hippoboscoidea: Streblidae and Nycteribiidae). Appl Environ Microbiol 2013; 79: 2952–2961.
104. Wilkinson DA, Duron O, Cordonin C, Gomard Y, Ramasindrazana B, Mavingui P, et al. The bacteriome of bat flies (Nycteribiidae) from the Malagasy Region: a community shaped by host ecology, bacterial transmission mode, and host-vector specificity. Appl Environ Microbiol 2016; 82: 1778–1788.
105. Graciolli G, Dick CW. Checklist of World Nycteribiidae (Diptera: Hippoboscoidea). https://www.researchgate.net/publication/322579074_CHECKLIST_OF_WORLD_NYCTERIBIIDAE_DIPTERA_HIPPOBOSCOIDEA. .
106. Graciolli G, Dick CW. Checklist of World Streblidae (Diptera: Hippoboscoidea). https://www.researchgate.net/publication/322578987_CHECKLIST_OF_WORLD_STREBLIDAE_DIPTERA_HIPPOBOSCOIDEA. .
107. Breitschwerdt EB, Kordick DL. Bartonella infection in animals: carriership, reservoir potential, pathogenicity, and zoonotic potential for human infection. Clin Microbiol Rev 2000; 13: 428–438.
108. Jiggins FM. Male-killing Wolbachia and mitochondrial DNA: selective sweeps, hybrid introgression and parasite population dynamics. Genetics 2003; 164: 5–12.
109. Hosokawa T, Koga R, Kikuchi Y, Meng X-Y, Fukatsu T. Wolbachia as a bacteriocyte-associated nutritional mutualist. Proc Natl Acad Sci U S A 2010; 107: 769–774.
110. Lack JB, Nichols RD, Wilson GM, Van Den Bussche RA. Genetic signature of reproductive manipulation in the phylogeography of the bat fly, Trichobius major. J Hered 2011; 102: 705–718.
111. Morse SF, Olival KJ, Kosoy M, Billeter S, Patterson BD, Dick CW, et al. Global distribution and genetic diversity of Bartonella in bat flies (Hippoboscoidea, Streblidae, Nycteribiidae). Infect Genet Evol 2012; 12: 1717–1723.
112. Nikoh N, Hosokawa T, Moriyama M, Oshima K, Hattori M, Fukatsu T. Evolutionary origin of insect-Wolbachia nutritional mutualism. Proc Natl Acad Sci U S A 2014; 111: 10257–10262.
113. Stuckey MJ, Chomel BB, de Fleurieu EC, Aguilar-Setién A, Boulouis H-J, Chang C-C. Bartonella, bats and bugs: A review. Comp Immunol Microbiol Infect Dis 2017; 55: 20–29.
114. Gibson CM, Hunter MS. Extraordinarily widespread and fantastically complex: comparative biology of endosymbiotic bacterial and fungal mutualists of insects. Ecol Lett 2010; 13: 223–234.