Pika Population Density Is Associated with the Composition and Diversity of Gut Microbiota.

Host population density is positively associated with the strength of social interactions or the frequency of physical contacts, and thus potentially influences microbial transmission among individuals. However, the relationship between host density and gut microbiota remains unknown. Here, we characterized the gut microbiota of plateau pikas (Ochotona curzoniae) in wild experimental fields with different host population densities. The abundance of some gut microbes significantly correlated with host density, such as Ruminococcaceae, Lachnospiraceae, and Staphylococcaceae. Intriguingly, host density was positively correlated with alpha diversity (Shannon diversity and observed species) of gut microbial communities. The inter-individual gut microbiota within high-density groups were more similar to each other than those of low-density groups. Host density significantly explained the variations in the gut microbiota, even after controlling sex, weight, diet and geographical locations. Based on the PICRUSt metagenome prediction, positive correlations were observed between host density and the relative abundances of 12 gene functions involved in cellular processes, environmental information processing and metabolism. These results indicate the importance of host density as a factor in shaping gut microbial composition and diversity in plateau pikas, and may further help us understand the social transmission of gut microbiota.

Introduction

Metacommunity theory in ecology is often used to describe the dynamics and composition of species assembly within a given patch (environment or habitat). In this framework, the neutral processes of species colonization, extinction and drift lead to the differentiation of communities across patches (Hubbell, 2001). Two patches that are close to each other may exhibit more similar communities due to the species colonization more frequently based on the higher probability of repeated species migration and re-colonization. When the metacommunity theory is applied to gut microbial communities, each host individual can be regarded as one patch. In this case, gut microbial communities are expected to be more similar among individuals in habitats with higher population densities, because the individuals from these habitats have more chances to come into contact with each other, and frequent contact may promote the horizontal transmission of microbial communities.

From the perspective of social evolution, social interactions can mediate exposure and susceptibility to bacteria, and socially mediated transmission is able to affect the evolutionary costs and benefits of social relationships (Lombardo, 2008; Archie and Theis, 2011). However, the influence of social relationships on the gut microbiota is an underappreciated outcome of group living, and is associated with disease transmission (Garrett et al., 2010; Henao-Mejia et al., 2012) and defense against pathogens (Koch and Schmid-Hempel, 2011). For instance, co-housing mice in the laboratory promotes the horizontal transmission of gut bacteria associated with obesity (Henao-Mejia et al., 2012) and inflammatory bowel disease (Garrett et al., 2010), indicating that the social interactions may increase the microbiota-associated disease risk. Moreover, socially transmitted bacteria protect bumblebees against a widespread and highly virulent gut parasite (Crithidia bombi) (Koch and Schmid-Hempel, 2011), suggesting that microbial transmission can provide some benefits for the host. It has been shown that host density was a determining factor of abundance in parasite communities, and that transmission rates of parasites rely on host population density (Arneberg et al., 1998). Similarly, the probability of microbial transmission through social interactions should depend closely on host population density in a limited space. A recent report showed that social interactions predict gut microbiome composition in wild baboons (Tung et al., 2015). However, measuring social interactions in the wild is quite difficult, potentially inaccurate and time-consuming. As an alternative, host population density is easy to measure and is associated with social interactions (Wang et al., 2005); thus it may potentially influence the composition of gut microbiota. Yet, to date, the relationship between host density and gut microbiota remains unknown.

Plateau pika (Ochotona curzoniae) is a keystone species for biodiversity and one of the most important animals shaping the landscape and function of grassland ecosystems in the Qinghai-Tibet Plateau (Smith and Foggin, 1999; Li et al., 2013). They are highly social herbivorous mammals that live in family groups (Dobson et al., 1998). Under moderate population density, they create microhabitats that result in an increase in plant species diversity and contribute positively to ecosystem-level dynamics (Smith and Foggin, 1999). Plateau pikas could reach a very high population density in the Qinghai-Tibet Plateau, especially in the degraded grassland. A high level of population density may also increase the transmission of gut microbes, and the risk of microbiota-associated diseases. However, it remains unknown whether pika population density is associated with gut microbiota. Here, we used plateau pikas as an animal model to explore whether host density is correlated with the composition and diversity of gut microbiota. Our results showed that host density is positively correlated with alpha diversity of gut microbiota, and the inter-individual microbiota within the high-density groups were more similar to each other than those within the low-density groups.

Materials and methods

Experimental design and sample collection

We chose two experimental fields (A and B) that harbored homogeneous grass environments and very similar plant communities. They were located in the near Zeku County of the Qinghai-Tibet Plateau. Due to the different extent of the killings by humans and the population migration of plateau pikas in the past few years, different population densities were naturally formed in these experimental sites. The geographical distance between A and B is about 80 kilometers. Each experimental field consisted of three 1 × 1 km plots. To reduce artificial disturbance, each plot was surrounded by a metal fence. Based on our long-term observations, there were almost no natural enemies of plateau pikas in these plots. The distance between the two plots in the same field was at least 500 m. The pika population density in each plot was measured based on the walked transect method (Qu et al., 2011). The population densities (individuals per hm2) of the six plots A1, A2, A3, B1, B2, and B3 were 55, 10.5, 2.5, 42.5, 13, and 3, respectively (Figure 1). In order to explore the possible differences of food resources for herbivorous pikas between fields A and B, the plant species were identified based on morphological characteristics, and the plant cover was measured in each experimental field according to the method described by Hogan (2010).

Figure 1

The experimental design in this study. A total of six plots from two experimental fields were selected. The geographical distance is 80 kilometers between two sites (A, B). The population densities of six plots, A1, A2, A3, B1, B2, and B3, were 55, 10.5, 2.5, 42.5, 13, and 3 individuals per hm2, respectively. All the plots can be divided into three groups, including high-density groups (A1 and B1), middle-density groups (A2 and B2), and low-density groups (A3 and B3).

All animal care procedures were consistent with the provision of the Institution of Animal Care and the Ethics Committee of Chengdu Institute of Biology, Chinese Academy of Sciences. Processing of wild animals and sample collection were strictly in accordance with the guidelines of our academic institutions. We randomly captured 7–9 individuals with ropes in each plot. After pikas were euthanized and dissected, their gut contents were immediately collected from the caecum of wild pikas, and stored immediately at −20°C in a portable refrigerator. A total of 50 caecum content samples were collected from the six plots, including A1 (n = 7), A2 (n = 8), A3 (n = 9), B1 (n = 9), B2 (n = 9), and B3 (n = 8). The corresponding information of each sample was listed in Table S1.

16S rRNA sequencing and bioinformatics analysis

The Ezup genomic DNA extraction kit for soil (Sangon Biotech, China) was used to extract total DNA of gut contents. DNA concentration and purity were determined using the Nanodrop 2000 Spectrophotometer. The extracted DNA was diluted to 10 ng/μL for PCR amplification. In order to amplify the V4–V5 hypervariable region of the microbial 16S rRNA gene, the universal primers 515F (5′-GTGYCAGCMGCCGCGGTA-3′) and 909R (5′-CCCCGYCAATTCMTTTRAGT-3′) with 12 nt unique barcode at 5′-end of 515F (Tamaki et al., 2011) were used. PCR amplifications were performed in duplicate with 25 μL reaction mix containing 1 × PCR buffer, 1.5 mM MgCl2, each deoxynucleoside triphosphate at 0.2 mM, each primer at 1.0 μM and 0.25 U of Ex Taq (TaKaRa, Dalian) and 10 ng genomic DNA. The thermal cycling procedure consisted of an initial denaturation step at 94°C for 3 min, followed by 30 cycles of 94°C for 40 s, 56°C for 60 s, and 72°C for 60 s, and a final extension at 72°C for 10 min. After PCR amplification, the two PCR products were mingled together and subjected to electrophoresis using 1.2% agarose gel. The correct band (~400 bp) was excised and purified using SanPrep DNA Gel Extraction Kit (Sangon Biotech, China) and quantified with the Nanodrop 2000 Spectrophotometer. All samples were pooled together with an equal molar amount from each sample. The sequencing samples were prepared using the TruSeq DNA kit according to the manufacturer's instruction. The purified library was diluted, denatured, re-diluted, mixed with PhiX (equal to 30% of final DNA amount), as described in the Illumina library preparation protocols, and then applied to an Illumina Miseq platform for sequencing (Reagent Kit V2) at the Environmental Genome Platform of Chengdu Institute of Biology.

The raw sequence data were processed using QIIME Pipeline-Version 1.7.0 (http://qiime.org/tutorials/tutorial.html). All sequences were trimmed and assigned to each sample based on their barcodes (barcode mismatches = 0). The overlapping paired-end reads were merged using the FLASH-1.2.8 software (Magoc and Salzberg, 2011). The merged sequences with high quality (reads length >300 bp, without ambiguous base “N,” and average base quality score >30) were used for further analysis. Due to possible contamination of chloroplast sequences in PCR amplification, we employed the Metaxa2 software tool to remove chloroplast sequences from our large sequencing datasets (Bengtsson-Palme et al., 2015). The aligned 16S rRNA gene sequences were used for a chimera check using the Uchime algorithm (Edgar et al., 2011). All sequences were clustered into operational taxonomic units (OTUs) at a 97% identity threshold using CD-HIT (Li and Godzik, 2006). Singleton OTUs were filtered out. To compare samples with different sequencing depths, each sample was rarefied to the same number of reads (8234 sequences). The alpha diversity indices, including observed species and Shannon diversity, were calculated. To assess beta diversity, the unweighted and weighted UniFrac distance metrics, which use phylogenetic information to calculate community similarity (Lozupone and Knight, 2005), were produced through the QIIME pipeline. In addition, in order to compare community similarity based on the presence/absence of OTUs, the Jaccard index was used; to compare community similarity based on OTU abundance, the Bray-Curtis dissimilarity matric was used with Hellinger transformed data (Dixon, 2003). The rarefaction curves were generated from the Goods coverage at the OTU level. Taxonomy was assigned using the Ribosomal Database Project classifier (Wang et al., 2007).

The original sequence data are available at the European Nucleotide Archive by accession number PRJEB11203 (http://www.ebi.ac.uk/ena/data/view/PRJEB11203).

Statistical analysis

Because the sample size was uneven in each group, an analysis of similarity (ANOSIM), which does not allow for interaction terms (Dill-McFarland et al., 2015), was applied to evaluate if microbial communities were significantly different across groups. ANOSIM analysis was implemented using the procedure “anosim” in the R “vegan” package (Warton et al., 2012). We calculated the average weighted UniFrac dissimilarity within each group. In order to explore the relationship between host density and the average weighted UniFrac dissimilarity, the Spearman correlation analysis was tested. Differences in the levels of variation within each group were tested with permutational tests of multivariate dispersions (PERMDISP). The Spearman correlations between the mean relative abundance (>0.05%) of microbial taxa at phylum, family and genus levels, and host population density were analyzed using SPSS 13.0 software. False discovery rate (FDR) values were evaluated using the Benjamini-Yekutieli method (Benjamini and Yekutieli, 2001).

Network analysis

In order to better understand the interactions among different species in the gut microbial communities and their correlations with the host population density, we constructed phylogenetic molecular ecological networks (pMENs) through Random Matrix Theory (RMT)-based methods, according to the molecular ecological network analyses pipeline (http://ieg2.ou.edu/MENA/main.cgi) with default parameters based on the OTUs from the 16S rRNA gene sequences. This approach is automatically defined and robust to noise, thus it provides remarkable solutions to several common issues associated with high-throughput metagenomic data (Deng et al., 2012). Only those OTUs identified in more than half of all individuals were used for network construction. The OTU modules were detected by fast greedy modularity optimization. The modules with more than 6 nodes were used for network-Eigengene analysis. The module-EigenGene analysis was then used to detect the correlations between microbial modules and host density. Cytoscape 3.0.2 software was used to visualize the final ecological network.

Predicted metagenomes

PICRUSTv1.0.0 (Langille et al., 2013) was used to predict the abundances of KOs from OTU abundances rarefied at 8234 reads per sample. We calculated the mean relative abundance of gene functions in each group. Thereafter, the spearman correlation between population density and predicted KEEG functions was analyzed using SPSS 13.0 software. P-values have been corrected using the false discovery rate control.

Results

Diversity of gut microbiota in plateau pikas

To compare samples with different sequencing depths, each sample was rarefied to 8234 sequences. At a threshold of 97% sequence identity, 37,258 unique OTUs were identified in all samples. The OTU-level rarefaction curve of Goods coverage across all samples has reached stable values (Figure S1), showing that most of the gut microbial diversity had already been captured in our study, despite the possibility of detecting rare new OTUs with an additional sequencing depth.

Plateau pikas harbored diverse gut microbial communities. Across all samples, 90.1 and 47.2% of total sequences were assigned to 29 phyla and 147 families, respectively. Firmicutes (mean ± SEM = 56.5% ± 1.4%), Bacteroidetes (30.2% ± 1.3%), Spirochaetes (1.3% ± 0.3%), and Proteobacteria (0.9 % ± 0.1%) were the dominant phyla across all samples. The composition of each sample was visualized at the phylum level in Figure 2. At the family level, the gut microbiota of the plateau pikas was dominated by Ruminococcaceae (18.1% ± 0.7%), S24-7 (13.4% ± 1.2%), Prevotellaceae (7.4% ± 0.5%) and Lachnospiraceae (2.8% ± 0.1%). However, only a small proportion of sequences (19.3 %) were identified at the genus level. The dominant genera were Prevotella (7.2% ± 0.5%), Oscillospira (5.1% ± 0.3%), Ruminococcus (2.1% ± 0.1%), and Treponema (1.1% ± 0.3%). The compositions of gut microbiota at phylum, family and genus levels were provided in Table S2.

Figure 2

The composition of pika gut microbiota at phylum level. Only those phyla with the relative abundances >0.1% were shown.

Population density is associated with the composition of pika gut microbiota

The Spearman correlation analysis was implemented to investigate the relationships between the relative abundances of major taxa (phyla, families, and genera) and host population density. At the phylum level, Acidobacteria were positively correlated with the host density, while Spirochaetes showed negative correlations with the host density. At the family level, Ruminococcaceae, Lachnospiraceae, and Staphylococcaceae were positively associated with population density, while Spirochaetaceae, Coriobacteriaceae, and Clostridiaceae showed negative correlations with the host density. At the genus level, Oscillospira, Roseburia, and Clostridium showed positive correlations with the host population density, whereas Treponema showed negative correlations with the host density (Table 1).

Table 1

Gut microbial taxa significantly correlated with pika population density.

	r	P-value
PHYLA
Spirochaetes	−0.943	0.005**
Acidobacteria	0.943	0.005**
FAMILIES
Ruminococcaceae	0.829	0.042*
Lachnospiraceae	0.886	0.019*
Spirochaetaceae	−0.943	0.005**
Coriobacteriaceae	−0.943	0.005**
Clostridiaceae	−0.886	0.019*
Staphylococcaceae	0.886	0.019*
GENERA
Oscillospira	0.829	0.042*
Treponema	−0.943	0.005**
Roseburia	0.829	0.042*
Clostridium	0.986	0.0003***

Spearman correlation analysis is used. P-value has been corrected using the false discovery rate control. The significance is indicated by

P < 0.05,

P < 0.01,

P < 0.001.

Population density is associated with the diversity of gut microbiota

The alpha diversity of pika gut microbiota from different groups was calculated (Table 2). Both the Shannon diversity (Spearman rho = 0.87, P < 0.05) and observed species (Spearman rho = 0.829, P < 0.05) showed significant positive associations with the host density.

Table 2

The alpha diversity of pika gut microbiota in different experimental plots and Spearman correlations with host population densities.

Experimental plots (population density: individuals/hm2)	Shannon diversity	Observed species
A1(55)	10.7 ± 0.1ab	3337 ± 103ab
A2(10.5)	10.6 ± 0.1ab	3330 ± 112ab
A3(2.5)	10.5 ± 0.13ab	3147 ± 91ab
B1(42.5)	10.8 ± 0.04a	3433 ± 23a
B2(13)	10.6 ± 0.1ab	3299 ± 73ab
B3(3)	10.4 ± 0.23b	3105 ± 204b
Spearman's correlation coefficients	0.87*	0.829*

Data are expressed as mean ± SEM. Different letters indicate significant differences of alpha diversity among experimental plots. The significant correlation is indicated by

P < 0.05.

In order to evaluate the beta diversity of pika gut microbiota, we calculated the Jaccard, Bray-Curtis, unweighted UniFrac and weighted UniFrac distance metrics across samples. The results using these distance metrics were very similar; thus, only the results for the weighted UniFrac were presented here. The gut microbiota in pikas were significantly influenced by the host density (ANOSIM, r = 0.2, P < 0.05). The microbial communities within the high-density groups (A3 and B3) or middle-density groups (A2 and B2) had less inter-individual variations than those of the low-density groups (A1 and B1). That is, the individuals within the high-density groups or middle-density groups clustered more tightly together than those within the low-density groups, Figure 3). The results were confirmed by PERMDISP (P < 0.001). In addition, we calculated the average weighted UniFrac dissimilarities of gut microbiota within each plot (Figure 4). The inter-individual gut microbiota within the high-density groups or middle-density groups were more similar to each other than those within the low-density groups (One-way ANOVA, P < 0.001). Host density was negatively correlated with the average inter-individual dissimilarity within each plot (Spearman rho = −0.886, P < 0.05). Notably, there was one outlier sample shown in Figure 3. However, we still obtained similar results even after removing this sample.

Figure 3

The Principal coordinate analysis of pika gut microbiota of high-density groups (A1 and B1), middle-density groups (A2 and B2), and low-density groups (A3 and B3).

Figure 4

The inter-individual variations of gut microbiota within each group. The significance is indicated by different letters (P < 0.001, one-way ANOVA).

Differences of gut microbial compositions among different population densities were unlikely explained by diet

Previous research showed that diet played an important role in shaping the differences of gut microbial compositions in primates and humans (Degnan et al., 2012; Yatsunenko et al., 2012). However, the two experimental fields in this study harbored a very homogeneous environment and similar plant communities. About 50% of plant coverage in the experimental plot A or B was contributed by Kobresia humilis, and similar proportions were contributed by other plant species, such as Ajuga lupulina, Leontopodium nanum, and Potentilla anserine (Figure S2). Only three plant species (Stipa capillata, Lagotis Gaertn, and Aconitum szechenyianum) were different between plots A and B (Table S3). Moreover, we calculated the effects of host density, diet diversity, body weight, geographical locations (A and B) and sex in shaping the gut microbiota of plateau pikas (Table S4). The results showed that only host density significantly influenced the gut microbial composition (ANOSIM, r = 0.2, P < 0.05), while diet diversity, host weight, sex and geographical locations did not show significant effects (ANOSIM, P > 0.05, Table S4).

Network analysis showed that density positive-correlation microbes were mainly Prevotella, Ruminococcus, and Oscillospira

Network analysis can reveal microbial interactions of different species in a microbial community, such as competition and mutualism. The species that own a similar niche or function form a module (Deng et al., 2012). Based on network analysis, the gut microbiota of plateau pikas consisted of 14 modules (Figure 5). Using module-EigenGene analysis, we found that the module 10, which mainly consisted of OTUs affiliated with Clostridiales, was negatively associated with host density, whereas modules 2, 3, 4, and 9, which mainly contained OTUs associated with Prevotella, Ruminococcus, and Oscillospira, showed positive correlations with host density (Table S5). The taxonomy of all the OTUs from each module was listed in Table S6.

Figure 5

The Molecular Ecological Networks (MEN) of pika gut microbial communities. Number in the circles is module number. Red line means positive correlation, and gray line means negative correlation.

Population density is associated with gene functions of gut microbiota

Correlations between host density and gene functions based on the PICRUSt metagenome prediction were tested. The relative abundances of 12 gene functions involved in cellular processes, environmental information processing and metabolism were positively correlated with host density (Table S7). In particular, those genes involved in glycan biosynthesis and metabolism, such as N-glycan biosynthesis, glycosaminoglycan degradation and other glycan degradation, were positively associated with host density (P < 0.05).

Discussion

Our results provided robust evidence that host population density is associated with the composition and diversity of gut microbiota in wild plateau pikas. To our knowledge, this is the first report to demonstrate that the host density can explain the variations in gut microbial communities after excluding the effects of sex, weight, diet, and geographical location. Host density was positively correlated with the strength of social interactions or the frequency of physical contacts. Thus, our results indicated that physical contacts or social relationships may play important roles in shaping gut microbial composition and diversity in natural populations.

Microbial transmission is very important in maintaining the diversity of the gut microbiota in wild mammals (Lombardo, 2008). Similar to mice and bumblebees (Garrett et al., 2010; Koch and Schmid-Hempel, 2011), plateau pikas are also coprophagic, thus a higher host density may increase the probability and incidences of microbial fecal-oral transfers among individuals. In addition to fecal-oral transfer, the proximity and physical contact between hosts was associated with the exchange of gut microbiota within the host species (Yatsunenko et al., 2012; VanderWaal et al., 2014). In experimental mice, hosts reared in the same cage are inclined to share their gut microbes (Campbell et al., 2012). In wild baboons, closer grooming partners harbor more similar gut microbial communities (Tung et al., 2015). For plateau pikas, the frequency of affiliative behaviors (e.g., kissing and contact) is closely correlated with host density (Wang et al., 2005), and potentially influences the transfer frequency from hand to mouth and from mouth to mouth. Social microbial transfer is not limited to the same animal species, but also occurs between different species. In addition to direct transmission, the microbial transmission can be indirectly mediated through shared environments, such as soil (Nunn et al., 2011).

Previous research on social-mediated transmission focuses on parasites (Nunn et al., 2011) or pathogens (Fournier et al., 2015), which may exert diseases or adverse effects on hosts. However, this study implicated that social transmission microbes may also provide powerful benefits for the host. In our study, the relative abundances of Ruminococcaceae and Lachnospiraceae were positively correlated with host density, and these density positive-correlation bacteria were likely social transmission microbes. Many members of these families harbor various genes encoding cellulase, hemicellulase or oligosaccharide-degrading enzymes (Dai et al., 2014); thus, they may play important roles in plant polysaccharide degradation. In particular, the genus Roseburia, one common butyrate producer (Duncan, 2002), showed positive correlations with host density, implicating that SCFA-producing bacteria can also be transmitted among pika populations. Intriguingly, the dominant genus Oscillospira, which was associated with fresh forage diets (Mackie et al., 2003), also showed a positive correlation with host density. Plateau pikas often consume a large amount of fresh grass in summer for their development and reproduction. When the host density increases in a limited space, the plateau pikas may compete for fresh grass more strongly, thus, increase the abundance of Oscillospira associated with fresh diets. In this study, metagenome prediction based on PICRUSt showed the relative abundances of function genes involved in glycan biosynthesis and metabolism were positively correlated with host density, indicating that host density may indirectly influence the function of gut microbiota. Future work may address the relationship between gut microbiota function and host population density in order to understand overall functional ecology of microbial transmission.

Some bacteria, such as Treponema, were negatively correlated with host density. However, these bacteria may have similar spatial and nutritional niches with those density positive-correlation microbes. Their abundances may be strongly influenced by the density positive-correlation microbes through microbial interactions. Understanding the function of these density negative-correlation microbes has important implications for the study of host health and nutrition. On the one hand, some members of the genus Treponema may cause mammalian diseases, although they are considered as nonpathogenic members in the gut community of ruminants and nonhuman primates (Stanton and Canale-Parola, 1979; McKenna et al., 2008). Future research is needed to explore whether these microbes are associated with the occurrence of pika disease. On the other hand, some species in this genus are associated with the digestion of specific soluble fibers (Bekele et al., 2011), and facilitate the cellulose degradation with co-occurring bacteria (Kudo et al., 1987).

The assumption that host density is associated with the alpha diversity of gut microbiota by positively affecting microbial transmission rates is crucial in epidemiological theory. We can infer that the transmission rates of microbes are accelerated from a minority of individuals to a wider range of host individuals in the high-density groups. As diversity is important for the stability and performance of all ecosystems (Clarke et al., 2014), the positive correlations between host density and alpha diversity of gut microbiota suggest that social microbial transmission may improve stability and performance of gut microbiota in high-density groups. Social transmission microbes in a high-density group may be more important to regulate inter-individual differences of gut microbiota. The hypothesis is consistent with our results that the inter-individual microbial compositions within high-density groups were more similar to each other when compared to those within low-density groups. Our findings supported that gut microbiota can be regarded as a result of “supply-side” ecology (Lewin, 1986); that is, you get what you are exposed to Bush (1990). Host density is an important determinant of exposure rates (Arneberg et al., 1998); thus higher microbial diversity and less inter-individual variations of gut microbiota in high-density groups are the results of higher exposure rates. From the view of the metacommunity theory, plateau pikas can be regarded as free-living animals occupying a set of discrete patches, and patch occupancy (i.e., prevalence of microbial species) often increases with patch density. It seems that low-density patches are less frequently occupied because they are less likely to come in contact with other pikas/microbes. Theoretically, the connection is formed by host density positively affecting microbial transmission rates. When host densities increase, those exposed bacteria (e.g., feces-associated bacteria or skin-associated bacteria) have more chances to contact with other host individuals.

Here, we highlighted the effects of transmission rates on gut microbiota, and considered contact rates as a sole correlation parameter of host density. However, we can't ignore the possible influences of other factors relevant to microbial transmission, such as host family size, kinship relationship, sociality and population structure. It remains uncertain whether population density has a widespread effect on microbial diversity in other host species. It would be particularly valuable for understanding the linkage between microbes-associated disease risk and host population density.

Author contributions

HL, JQ designed the experiments. HL, JQ, TL, and QL performed the field experiments. HL, XL wrote the paper. JL revised the paper.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.