The sampling site, Xianren Cave (102°20′E, 24°30′N), is located on the Yunnan-Guizhou Plateau at an altitude of 2049 m in Yunnan Province, China. The site has a subtropical climate with an average annual temperature of 15.5°C and an average annual rainfall of 973.5 mm, with precipitation mainly concentrated in June–October. The vegetation surrounding the cave mainly included Juglans catanensis, Pinus yunnanensis furnace, Cryptomeria fortunei, and Phyllostachys sulphurea, as well as agricultural fields and villages, all of which provided a good foraging habitat for bats. Preliminary monitoring revealed that five rhinolophid bat species coexisted stably in the cave, including R. macrotis (about 50 individuals), R. osgoodi (about 40 individuals), R. ferrumequinum (about 20 individuals), R. affinis (about 20 individuals), and R. pusillus (about 20 individuals) (unpublished data).

To minimize the disturbance to lactating bats, we chose to conduct the study in August–September 2017 and 2018, from 21:00 to 4:00 the following morning. Bats were captured with mist nets at the entrance of the Xianren cave after they returned to the cave. Species and sex were identified based on external morphological characteristics, and whether bats were adults or juveniles was judged based on epiphyseal fusions (Stebbings, 1968). Once juveniles were identified, they were immediately released in situ. Captured bats were placed singly in sterilized kraft paper bags for 2 h while droppings were collected. The bags were checked frequently to ensure the freshness of the droppings. The droppings were placed into lyophilized tubes containing RNAlater (RNA-EZ Reagents RNA-Be-Locked A, Sangon Biotech, Shanghai, China) using sterilized forceps. At least five droppings per bat were collected, returned to the laboratory, and stored at −80°C for DNA extraction.

After the guano was collected, the forearm length of each bat was measured using a digital caliper (TESA-CAL IP67, Tesa Technology, Renens, Switzerland, 0.01 mm) and body mass was determined with an electronic balance (ProScale LC-50, Accurate Technology, Inc., Asheville, NC, United States, 0.01 g). The wing morphology of each bat was drawn with reference to Thabah’s method (Thabah et al., 2007). The bite force of the molars was recorded five times for each bat following the method of Freeman and Lemen (2008) using a bite force measurement system (Nanjing Bioinspired Intelligent Technology Co., Ltd., Nanjing, China) with an accuracy of 0.01 N. The system includes a sensor (NBIT-S1-100NHL-001) and data acquisition system (NBIT-Dus-2404A), and was used following the manufacturer’s protocol. The maximum bite force was taken as the bat’s bite force value. The echolocation calls were recorded using an ultrasound recorder (UltraSound Gate 116, Avisoft Bioacoustics, Berlin, Germany) when the bats were resting in a temporary laboratory (5 m × 4 m × 3 m), with the recording microphone aimed at the resting bat’s head at a distance of 30 cm. The sampling frequency was 375 kHz, and the resolution was 16 bits. The data were input into a laptop computer and analyzed using Avisoft SASLab Pro software. Each bat was recorded for at least 1 min (Jacobs et al., 2007; Li et al., 2007). After all measurements and acoustic recordings were finished, the bats were released in situ. A total of 98 bats were captured in the field, including 36 R. macrotis (18 females and 18 males), 29 R. osgoodi (21 females and eight males), 11 R. ferrumequinum (six females and five males), 12 R. affinis (six females and six males), and 10 R. pusillus (two females and eight males).

Genomic DNA was extracted from 200-mg guano per bat using the E.Z.N.A™ Mag-Bind Soil DNA Kit (OMEGA Bio-Tek, Norcross, GA, United States) according to the manufacturer’s instructions. After DNA extraction, DNA integrity was tested using 2% agarose gels. A 225-bp fragment of the cytochrome c oxidase subunit I (COI) was amplified using the primers LCO-1490 (5′-GGTCAACAAATCATAAAGATATTGG-3′) and ZBJ-ArtR2c (5′-WACTAATCAATTWCCAAATCCTCC-3′) (Zeale et al., 2011; Chang et al., 2019). Genomic DNA was accurately quantified using the Qubit 2.0 DNA Assay Kit (Life Technologies, Thermo Fisher Scientific, Waltham, MA, United States) to determine the amount of DNA added to the PCR reaction. The first step of PCR amplification reaction was conducted in a final volume of 30 μl containing 15 μl 2 × Taq Master Mix, 10 ng of genomic DNA, 1 μl primer F, and 1 μl primer R, with water added to reach a volume of 30 μl. The conditions for PCR were as follows: 94°C for 3 min; five cycles at 94°C for 30 s, 45°C for 20 s, and 65°C for 30 s; 20 cycles at 94°C for 20 s, 55°C for 20 s, and 72°C for 30 s; and a final extension at 72°C for 5 min. After the first step, the PCR products were checked using 2% agarose gel. The second step of the PCR amplification reaction was conducted in a final volume of 30 μl containing 15 μl 2 × Taq Master Mix, 10 ng of genomic DNA (from the first step), 1 μl primer F, 1 μl primer R, and water added to reach a volume of 30 μl. The PCR amplification reaction conditions were as follows: 95°C for 3 min; five cycles at 94°C for 20 s, 55°C for 20 s, and 72°C for 30 s; and a final extension at 72°C for 5 min. The PCR products were purified using Agencourt AMPure XP (Beckman Coulter, United States), and the DNA concentration of each sample quantified using the Qubit 2.0 DNA Assay Kit (Life Technologies, Thermo Fisher Scientific, Waltham, MA, United States) to normalize samples according to the manufacturer’s protocol. The final products were sequenced using the Illumina Miseq platform (2 × 300 bp) at Sangon Biotech (Shanghai, China).

Primers and adaptors were removed from the raw sequences using cutadapt v1.2.1 (Martin, 2011), and paired-end reads were merged using PEAR v0.9.6 (Zhang et al., 2014) based on barcode tags to distinguish samples. Finally, the data files were quality-filtered using Prinseq v0.20.4 (Schmieder and Edwards, 2011). The operational taxonomic units (OTUs) were clustered using Usearch after the singletons and chimeras were removed (Edgar, 2010). All optimized sequences were mapped to representative sequences, and OTU tables with a 97% similarity threshold were generated. The samples with fewer than 12,000 sequences were discarded, and the data were rarefied to 12,000 sequences per sample. Following this standard, six individuals with fewer than 12,000 sequences were removed (three R. macrotis, one R. osgoodi, one R. ferrumequinum, and one R. affinis). Finally, to minimize the effect of sequencing errors, the OTUs representing <0.1% of the normalized sequences for each sample were removed to prevent the generation of potentially erroneous results (Bokulich et al., 2013). Representative sequences of each OTU were compared with the reference sequences in the Barcode of Life Database (BOLD¹) and the Genbank database² to obtain taxonomic information. OTUs that did not match any taxonomic information were excluded. The identification criteria were based on slightly modified “strict” and “best” matching methods (Ross et al., 2008; Razgour et al., 2011).

After scanning the mapped wing shape, the wing span and wing area of each individual were measured using AutoCAD 2021 (Autodesk, United States), and the aspect ratio and wing loading were calculated for each individual (Norberg and Rayner, 1987; Salsamendi et al., 2012). For each individual bat, after removing the first and last syllables of each echolocation call sequence, 15 syllables with a signal-to-noise ratio higher than 40 dB were randomly selected. Due to the special cochlear structure and Doppler-shift compensation of rhinolophid bats during flight, bats are able to actively adjust their call frequency to ensure the dominant frequency keep stable (Russell and Kossl, 1999; Metzner et al., 2002; Davies et al., 2013). Therefore, three acoustic parameters, including peak frequency (PF), minimum frequency (Fmin), and maximum frequency (Fmax) were analyzed on the power spectrogram and oscillogram with Avisoft-SASLab Pro version 5.2.9 (Avisoft Bioacoustics, Berlin, Germany) based on a 1,024 FFT, 100% frame size, and 93.75% temporal overlap.

We quantified the diet composition using the percent of occurrence (POO) and relative read abundance (RRA) at the order level (Deagle et al., 2018). Chi-square tests (SPSS 22.0) and Kruskal–Wallis tests (R 4.1.1) with a post hoc Dunn test were used to test for differences in the prey orders consumed between bat species. The dietary niche breadth of each bat species was calculated using a standardized Levins’ index (B_A) (Razgour et al., 2011), the diet composition overlap between bat species was measured using Pianka’s measure of niche overlap (O_jk) in the “EcoSimR” package in R (Gotelli et al., 2015), and the “qgraph” package was used to generate a network graph of dietary overlap between species. The ordination of samples was visualized depending on dietary composition at the OTU level with non-metric multidimensional scaling (NMDS), and an analysis of similarities (ANOSIM) test was used to determine whether the Bray–Curtis distance of dietary composition at the OTU level was greater between than within bat species (Oksanen et al., 2020). The ANOSIM test generates a global R statistic that provides a measure of the magnitude of difference among groups, with a range between −1 and + 1. When R > 0, the difference between groups was greater than the difference within each group. If the test results were significantly different overall, then a pairwise analysis between species was performed, and the p-values were corrected with Bonferroni correction. The differences between species were calculated for forearm length, body weight, aspect ratio, wing loading, bite force, and echolocation call parameters using one-way ANOVA or non-parametric Kruskal–Wallis tests depending on the test of normal distributions or the homogeneity of variances according to Shapiro and Bartlett tests in the package stats. The p-values for multiple comparisons were corrected using Bonferroni correction. All the data are presented as the means ± SD. Spearman correlation analysis was conducted to test the relationship between bite force and weight. To determine the difference in the alpha diversity of the diet among five rhinolophid bats at the OTU level, Kruskal–Wallis tests were used to compare the values of the Shannon–Wiener diversity and Gini–Simpson diversity with non-normal distributions in the “vegan” package. The effects of body size, wing morphology, bite force, and echolocation call parameters on the diet composition of bats were investigated using Canonical correspondence analysis (CCA) with the “vegan” package. The response variable (dietary species abundance data) and explanatory variables (body size, wing morphology, bite force, and echolocation calls) were Hellinger transformed and logarithmically transformed, respectively. The variance inflation factors (VIF) were calculated to identify the collinearity between the explanatory variables until the values of all factors were less than 10. To distinguish the contribution of a single variable, the relative importance of each influencing factor independently accounting for the variations in diet composition was quantified by applying a hierarchical partitioning analysis using the “rdacca.hp” package (Lai et al., 2022). All statistical analyzes and visualizations were conducted in R 4.1.1 (R Core Team, 2021).

A total of 7.64 million raw sequences were obtained, with an average of 83,061 sequences per sample. Following quality control, 6.19 million sequences were obtained, with an average of 67,303 sequences per sample. After deleting OTUs that did not match the taxonomic information, 1,031 OTUs were obtained through comparison with the BOLD and Genbank databases, all of which were from the Arthropoda phylum. Among all OTUs, 83.41% were identified at the order level, 45.59% at the family level, 19.79% at the genus level, and 14.06% at the species level (Supplementary Table S1). Insects comprised the majority (98.81%) of phylum detected. Lepidoptera (58.65%) and Diptera (21.72%) comprised the majority of all 15 orders detected. Noctuidae (20.21%), Geometridae (14.04%), Erebidae (10.00%) and Tortricidae (4.68%) comprised the majority of families detected. Some differences were detected in the diet composition among bat species. The diets of R. macrotis and R. ferrumequinum were mainly composed of Lepidoptera (R. macrotis: 96.31%, R. ferrumequinum: 89.88%). The diet of R. osgoodi consisted mainly of Lepidoptera (73.49%) and Diptera (17.28%). The diet of R. affinis mainly consisted of Diptera (59.39%) and Lepidoptera (26.64%). The diet of R. pusillus mainly consisted of Diptera (75.13%) and Coleoptera (15.40%).

At the OTU level, no significant differences were found in the Shannon–Wiener diversity or Gini–Simpson diversity between five rhinolophid bats (p = 0.94 and p = 0.91, respectively; Supplementary Figure S1). Overall, low levels of dietary overlap were found among bat species, with the greatest overlap value between R. affinis and R. pusillus (0.33), and the lowest value between R. macrotis and R. affinis (0.01; Table 1; Figure 1). The Levins’ standardized measure of niche breadth showed that dietary niche breadth ranged from 0.42 to 0.76, with the largest value for R. ferrumequinum and the smallest for R. macrotis (Table 1). The RRA analysis of the diet composition of each bat species at the order level showed significant differences in the consumption of Coleoptera, Diptera, Lepidoptera, Hemiptera, and Neuroptera (Kruskal–Wallis test, p < 0.05; Table 2; Figure 2). In addition, analysis of POO at the order level showed significant differences in the consumption of Coleoptera, Diptera, Hemiptera, and Neuroptera (Chi-squared test, p < 0.05; Table 2; Figure 2).

Body size differed greatly among bat species (Table 3; Supplementary Table S2), with the largest body weight and forearm length found in R. ferrumequinum (body weight: 18.58 ± 0.97 g; forearm length: 61.17 ± 0.40 mm) and the smallest found in R. pusillus (body weight: 4.90 ± 0.34 g; forearm length: 39.18 ± 1.21 mm). The maximum bite force was also found in R. ferrumequinum (5.65 ± 0.51 N), and the minimum was found in R. pusillus (1.02 ± 0.14 N). A significant positive correlation was detected between bite force and the body weight of bats (Spearman r = 0.87, p < 0.001). After removing the effect of body weight, there was still a significant difference between the four paired groups (R. ferrumequinum–R. macrotis, R. ferrumequinum–R. osgoodi, R. ferrumequinum–R. pusillus, and R. affinus–R. osgoodi) (Kruskal–Wallis test, p < 0.05; Table 3; Supplementary Table S2). R. macrotis had the lowest peak frequency (57.88 ± 0.86 kHz) and R. pusillus had the highest peak frequency (101.92 ± 1.49 kHz; Table 3; Figure 3). Significant differences were detected between bat species in body weight, forearm length, wing loading, bite force, peak frequency, minimum frequency, and maximum frequency (Kruskal–Wallis test, p < 0.001), as well as in the aspect ratio among species (one-way ANOVA, F = 5.16, p < 0.001; Table 3; Supplementary Table S2).

The NMDS ordination of diet composition between groups of five rhinolophid bats resulted in a two-dimensional solution with final stress of 0.103. The samples obtained from the five rhinolophid bats were spread out in the diagram, and the ANOSIM analysis showed significant differences in dietary composition in all pairwise comparisons (Figure 4A; r = 0.17, p = 0.001), except between R. affinis and R. pusillus (r = 0.09, p = 0.07). Because all factors had VIF values of less than 10 (Table 4), all factors were added to the model. The results of the CCA showed a significant relationship between all factors and diet data (adjusted r² = 1.63%, p = 0.001; Table 4; Figure 4B), and the results of the permutations tests for each influencing factor individually showed that body weight, forearm length, peak frequency were all significantly correlated with dietary composition (p < 0.05; Table 4; Figure 4B). Hierarchical partitioning analysis showed that echolocation calls (peak frequency) were the most important factors, followed by body size factors (forearm length and body weight), while parameters such as wing loading and bite force had no significant effect on the diet composition structure for each bat species (Table 4; Figure 4B).

All five rhinolophid bats investigated in this study relied on Lepidoptera as their main diet type. This finding was similar to the conclusions of previous studies conducted in summer [R. ferrumequinum (Jin et al., 2005; Andreas et al., 2013); R. affinis (Ye et al., 2009); R. pusillus (Wei et al., 2006); R. osgoodi and R. macrotis (Shi et al., 2009)]. However, there are some differences between previous findings and the current study. For example, studies on R. affinis and R. pusillus showed that Coleoptera were an important diet component (R. affinis: 39.51%, R. pusillus: 41.08%; Wei et al., 2006; Ye et al., 2009), while the present study found that Coleoptera make up only a relatively small part of the diet (R. affinis: 5.14%, R. pusillus: 15.40%), which may be due to the difference in available insect resources at different study sites or different research methods used. These insects provide a substantial, high-quality food source for bats that can ensure that their nutritional needs are met (Clare et al., 2011; Emrich et al., 2014; Rolfe et al., 2014; Chang et al., 2019). The present study used high-throughput sequencing to identify the species of insect residues remaining in bat guano, yielding a total of 145 identified insect species. This method has a high resolution compared with traditional guano analysis methods. For instance, this study detected small, soft-bodied Chironomidae in bat guano, which are easily overlooked in traditional guano analysis methods (Krüger et al., 2014).

Trophic niche differentiation is an effective way for sympatric bats to reduce competition for food resources and thus promote coexistence (Alley, 1982; Letten et al., 2017), and it has been demonstrated in numerous studies, such as Pteronotus macleayii and Mormoops blainvillii in Jamaica, and M. blainvillei, Pteronotus quadridens, and Pteronotus parnellii portoricensis on the island of Puerto Rico (Rolfe and Kurta, 2012; Emrich et al., 2014). In this study, all five rhinolophid bats had wide trophic niche breadth (0.42–0.76) and low overlap in trophic niche among species (0.01–0.33). This may have been due to the fact that sampling was conducted during the rainy season, when insects were abundant in our studied region (Zhang et al., 2017). However, significant differences on food composition among species (Figure 4A) could also suggest that morphological and acoustic differences shape their trophic relationships. In addition, according to optimal foraging theory, bats can prey upon insects based on their phenotypic characteristics; for example, bats with different peak frequencies will prey upon insects with different body sizes (Pyke et al., 1977; Perry and Pianka, 1997; Mello et al., 2004; Rollinson et al., 2013). In fact, significant differences were found in the insects captured by bats at the order level, such as Lepidoptera, which were the most abundant in the diet compositions of R. macrotis, R. osgoodi, and R. ferrumequinum. This may be because Lepidoptera are frequently a primary food resource of insectivorous bats due to their high abundance and wide global distribution (Hebert et al., 2003). This finding may also be due to the large size of lepidopterans, which have soft bodies and contain more energy, making them preferred prey for many predators (Krüger et al., 2014; Wray et al., 2021).

It has been shown that the trophic niche differentiation and stable coexistence of sympatric bats are influenced by many factors. In this study, differences were detected in the echolocation calls, body size, wing morphology, and bite force among five rhinolophid bat species. These differences may affect the flight patterns, foraging areas, and prey sizes of bats, and ultimately affect bat diet differentiation (Schoener, 1971; Norberg and Rayner, 1987; Mancina et al., 2012; Emrich et al., 2014).

Bats used echolocation calls to locate insects during foraging, and the parameters of echolocation calls, such as peak frequency, can significantly affect bat foraging habitats and prey types (Fullard et al., 1991; Bogdanowicz et al., 1997; Andreas et al., 2012; Gordon et al., 2019). In the present study, the results revealed that the peak frequency is a significant variance of all explanatory variables. This suggests that the echolocation calls of different bat species have an important influence on their dietary habits, similar to the findings of a previous study (Andreas et al., 2012). In general, bats with higher peak frequencies (short wavelengths) have smaller body size and fly at slower speeds, although they are more maneuverable than larger bats. These smaller bats are well-suited for foraging in dense spaces, and high-frequency echolocation calls help the bats move quickly and locate insects in cluttered environments (Jacobs et al., 2007). The prey detection hypothesis suggests that differences in echolocation frequencies may cause differentiation in prey sizes and that high-frequency calls are more sensitive to detecting smaller prey items (Barclay and Brigham, 1991; Jones, 1997). For example, R. euryale and R. mehelyi often exhibit sympatric distribution and are acoustically similar (R. euryale: 104.4 kHz; R. mehelyi: 106.8 kHz), but R. euryale tends to consume larger moths, while R. mehelyi prefers smaller moths (Salsamendi et al., 2006; Arrizabalaga-Escudero et al., 2018). A study on two cryptic species of pipistrelle bats found that 10 kHz difference in echolocation call frequency was not sufficient to influence the target strengths of the main prey types (Jones and Barlow, 2001). However, the large range of peak frequencies (57.88–101.92 kHz) detected among the bats in this study may have contributed to the differences in diet composition. For instance, in a previous study (Shi et al., 2009) conducted in the same cave as our study, R. macrotis had lower peak frequencies (~58 kHz) and a higher consumption of larger insects than coexisting R. lepidus (later identified as R. osgoodi, ~91 kHz) based on the measurement of prey size (Shi et al., 2009).

The morphological characteristics of bats are an important factor influencing the trophic niche differentiation of sympatric species. To some extent, morphological differences can reflect functional differences, which in turn affect the trophic niches of bat species (Pyke et al., 1977; Norberg and Rayner, 1987; Barclay and Brigham, 1991; Gnocchi et al., 2019). In the present study, body weight and forearm length were greater in R. ferrumequinum and R. affinis compared with the other species, with no significant difference between these two species and a significant difference from the remaining three species. The result of hierarchical partitioning analysis implies that forearm length and body weight play important roles in determining the trophic niche of sympatric bats. For instance, R. ferrumequinum and R. pusillus with different body sizes have significant differences in predation on Dipteran and Lepidopteran insects. Salsamendi et al. (2012) found that the niche partitioning of two sympatric sibling rhinolophid bats, R. euryale and R. mehelyi, was due to their differences in wing morphology, although their echolocation calls are almost identical. Andreas et al. (2013) suggested that the body size of three occurring horseshoe bats significantly affected the size of the prey. This suggests that, in addition to differences in echolocation calls, the body size of coexisting bats plays an important role in the differentiation of trophic niches, which is consistent with the predictions of morphological ecology (Findley and Black, 1983; Jacobs and Barclay, 2009).

Wing loading is an important parameter used to measure the flight ability of bats, and species with different wing morphology adapt to the local environment and food resources, allowing the more efficient use of food and habitat resources (Norberg and Rayner, 1987). In general, bats with high wing loading fly fast and forage in more open environments (Ye et al., 2009; Andreas et al., 2012; Mancina et al., 2012; Salsamendi et al., 2012). In this study, wing loading explained 15.0% of the diet composition, implying that wing loading had influence on inter-species trophic niche differentiation in bats. Not all pairwise comparisons of wing loading among bat species have shown significant differences, which may dilute the role of wing loading on interspecific trophic niche differentiation in bats. Furthermore, the bite force of bats is closely related to prey type, is influenced by factors including bat body weight and skull morphology, and is an important factor influencing the differentiation of trophic ecological of bats (Freeman and Lemen, 2010; García-Herrera et al., 2021). The bite force of bats in this study only explained a small part of diet composition, but was significantly correlated with body weight. This was similar to previous studies, in which larger bat species were generally found to have strong jaw and bite muscles that could produce greater bite force (Aguirre et al., 2003; Herrel et al., 2005; Freeman and Lemen, 2008). However, R. pusillus in this study consumed more Coleoptera than other species (RRA: 0.15; POO: 0.22) with small body size and low bite force, implying that bite force was not the main factor determining differences in the trophic ecology of bats. Because of its particularity (flight, nocturnal activity, migration, and social life), the trophic niche of a given bat species is dynamically affected by many factors, such as the spatial and temporal variation of food resources and environmental condition, as well as body condition, so the main factors affecting diet composition may change. It is difficult for a single factor to determine the trophic niche differentiation of bat species, and studying single factors alone may negatively impact our comprehensive understanding of the mechanisms of species coexistence in communities (Vesterinen et al., 2018; Gordon et al., 2019; Novella-Fernandez et al., 2020; Wray et al., 2021).

The CCA showed a low total explanation in this study. This pattern is more typical in ecological data (<10% of the variance explained) (ter Braak and Verdonschot, 1995) due to the nature of the presence–absence data type and the data including a large number of zeros, as was the case in this study. The dropping sample obtained from a single individual would not contain the full range of prey at the order level, but only one or a few arthropod orders. Although second-generation sequencing methods were used to obtain OTUs for species identification, this process generated more zeros in the data than traditional guano analysis methods.

In summary, this study investigated the trophic ecology and the effects of body size, wing morphology, bite force, and echolocation calls on the diet composition of five sympatric rhinolophid bats. The results revealed that all bats occupied a wide breadth of trophic niches in summer and that there was a low degree of trophic niche overlap between species. Significant differences in diet composition between species (except for between R. affinis and R. pusillus) were probably caused by the combined effect of the factors included in the study. CCA showed that body size and echolocation calls are the main factors that influenced the diet composition of bat communities, causing differentiation in trophic ecological niches among species, suggesting that the combined effects of multiple factors influenced the stability of bat communities.