Phylogenetic Community Structure and Niche Differentiation in Termites of the Tropical Dry Forests of Colombia

The mechanisms that structure species communities are still debated. We addressed this question for termite assemblages from tropical dry forests in Colombia. These forests are endangered and poorly understood ecosystems and termites are important ecosystem engineers in the tropics. Using biodiversity and environmental data, combined with phylogenetic community analyses, trait mapping, and stable isotopes studies, we investigated the termite community composition of three protected dry forests in Colombia. Our data suggest that the structuring mechanisms differed between sites. Phylogenetic overdispersion of termite assemblages correlated with decreasing rainfall and elevation and increasing temperature. Food niche traits—classified as feeding groups and quantified by δ15N‰ and δ13C‰ isotope signatures—were phylogenetically conserved. Hence, the overdispersion pattern implies increasing interspecific competition with decreasing drier and warmer conditions, which is also supported by fewer species occurring at the driest site. Our results are in line with a hypothesis that decreased biomass production limits resource availability for termites, which leads to competition. Along with this comes a diet shift: termites from drier plots had higher δ13C signatures, reflecting higher δ13C values in the litter and more C4 plants. Our study shows how a phylogenetic community approach combined with trait analyses can contribute to gaining the first insights into mechanisms structuring whole termite assemblages.


Introduction
The drivers that structure communities are still debated [1][2][3][4]. Neutral and/or deterministic mechanisms have been proposed to explain community assembly [5,6]. Neutral models highlight that mainly stochastic processes drive local communities. Species are regarded as ecologically equivalent. At larger scales, meta-communities are influenced by dispersal, speciation, and extinction [7]. On the other side, deterministic models describe local communities as an "arranged" assembly of species, based on their physiology and their defined niches [8,9]. These are two extreme views of processes affecting community composition. Real communities will often fall along a continuum containing components of randomness, as well as determinism. Whether processes differ systematically between taxa, habitats, geographic regions, or biomes is still unclear. A major unsolved question that remains is whether tropical ecosystems differ systematically from, for example, temperate regions, and whether such differences in the structuring mechanisms can contribute to explaining their high species richness. However, such species-rich ecosystems are notoriously difficult to study due to their high number of undescribed species. Two approaches can contribute to overcoming this hurdle. First, molecular barcoding uses short genetic markers in an organism's DNA to identify it as belonging to a particular species [10]. Hence, species can be identified more easily. Second, phylogenetic community analyses deduce potential mechanisms from the co-occurrences of species and their phylogenetic relatedness. They combine phylogenetic data with distributional and ecological data to assess whether and how communities of species differ from random assemblages with regard to evolutionary relatedness [1][2][3][4][5][6]. For instance, if species which locally co-exist are less closely related on the phylogeny than a random selection of species that could potentially co-occur (i.e., species from the regional species pool), this can indicate that interspecific competition can play an important role in structuring communities, given closely related species share the same niche traits [5,11]. As genetic sequence information and phylogenies are becoming increasingly available for many taxa, the use of phylogenetic community analyses can be a helpful and easy tool to gain first insights into potential assembly mechanisms. They also allow standardized comparisons between sites within a habitat, and between habitats, regions, or disturbed and natural areas, to inform about changes and similarities in community structuring mechanisms (for termites: [12,13]). Thus, both genetic barcoding and phylogenetic community analyses can help to gain first insights into community structuring mechanisms.
Tropical dry forests are the most threatened of all major tropical forest types [14,15]. Colombia has one of the best-conserved areas, mainly along the Caribbean coast [16]. These poorly studied ecosystems are threatened by land use, climate change [14], and urban expansion [15]. Termites (Termitoidea) are important ecosystem engineers of such tropical ecosystems [17,18]. They are important food sources for a wide range of species [19][20][21][22][23]. As the main macro-detritivores, they essentially contribute to the biotransformation of wood and litter into organic matter and the re-distribution of structural soil components [24]. Tropical forests produce plenty of dead plant material, which termites consume [25,26]. Four functional feeding groups are distinguished in termites [27,28]: dead wood-feeders (group I); dead wood, leaf, plant-litter feeders (group II); humus feeders (group III); and true soil feeders (group IV). No fungus-growing termites occur in the Neotropics, hence the differentiation of the feeding group IIF (i.e., fungus-feeder) is irrelevant [29]. Most termites of dry forests feed on twigs and litter and belong to feeding group I and II, while soil-feeders (sensu Anoplotermes-group) are relatively scarce in richness and abundance [26]. Nitrogen and carbon stable isotope ratios have been used in termites to elucidate feeding habits in more detail, including dietary preferences [30,31] and niche food differentiation [32][33][34][35]. However, most studies on stable isotopes have been conducted in savannas [31,33,35] and rainforests [32,36,37]. Information related to the trophic ecology of termites in Neotropical dry forests is unknown. We used isotope analyses to characterize the feeding niche of termites and combined this approach with phylogenetic community analyses to gain first insights into the mechanisms that may structure termite assemblages in dry tropical forests.
In a former study, we characterized the termite communities of these three sites by determining species diversity and abundances and associating them with environmental variables [26]. We studied fives transect belts per site (hereafter called study plots) using the standardized belt transect sampling protocols of Jones and Eggleton [43] and Hausberger and Korb [44] developed for termites. We surveyed each site by sampling a transect measuring 2 m × 100 m, divided into twenty 2 m × 5 m sections. Each section was searched for termites on the ground, and in trees, mounds, and soil, (eight soil pits 15 cm ×15 cm ×15 cm depth) for 30 min by two trained persons. All study plots were randomly chosen and they were separated from each other by, on average, around 560 m (min: 225, max: 1043, SD +/− 253 m) in Colosó, 1074 m (min: 366, max: 1982, SD +/− 557 m) in Ceibal, and 1606 m (min: 508, max: 3157, SD +/− 985 m) in Tayrona. We also took soil and litter samples and retrieved climate data from WorldClim v 1.4 (http://www.worldclim.org/). The data layers were generated through the interpolation of average monthly climate data from weather stations on a 30 arc-second resolution grid (often referred to as a "1 km 2 " resolution). Variables included were monthly total precipitation and monthly mean temperature (for more details see http://www.worldclim.org). A combination of morphological and genetic analyses (molecular barcoding) revealed a total of 32 species for all three sites (Table A1). . The data layers were generated through the interpolation of average monthly climate data from weather stations on a 30 arc-second resolution grid (often referred to as a "1 km 2 " resolution). Variables included were monthly total precipitation and monthly mean temperature (for more details see http://www.worldclim.org). A combination of morphological and genetic analyses (molecular barcoding) revealed a total of 32 species for all three sites (Table A1).

Determination of Food Niche
To determine the feeding type and characterize the food niche using δ 13 C and δ 15 N isotope analyses, we used specimens and material from a former study [26]. δ 13 C and δ 15 N isotope analyses were done for termite workers, but also for soil and leaf litter (leaf and small pieces of wood), which is potential food for the termites. For each termite species, a whole termite was used. As in several other studies [30,31,34,45,46], we could not exclude the gut as this would have left too small amounts to conduct the analyses. Prior to analyses, termite samples were stored in ethanol (>99.5% Merck, Darmstadt, Germany). Three replicates (if available) per site were analyzed. Only workers were taken into account to eliminate the effect of inter-caste differences in isotopic values, which could bias crossspecies comparisons [35,36]. In addition, workers are the caste that does the foraging and feeding within colonies. Five replicates of soil samples were collected from the top horizon (0-15 cm) at each site-one at each study plot-resulting in 15 samples in total for all study sites. Soil samples were cooled and directly dried after each field trip, and then sealed in plastic bags. Additionally, litter samples were collected on the ground; three samples were taken per study plot, including one at the start, in the middle, and at the end of each study plot, resulting in 15 replicates per site (45 in total). They included leaves, twigs, and dead wood. Like the soil samples, they were dried and kept cool, prior to analysis.

Determination of Food Niche
To determine the feeding type and characterize the food niche using δ 13 C and δ 15 N isotope analyses, we used specimens and material from a former study [26]. δ 13 C and δ 15 N isotope analyses were done for termite workers, but also for soil and leaf litter (leaf and small pieces of wood), which is potential food for the termites. For each termite species, a whole termite was used. As in several other studies [30,31,34,45,46], we could not exclude the gut as this would have left too small amounts to conduct the analyses. Prior to analyses, termite samples were stored in ethanol (>99.5% Merck, Darmstadt, Germany). Three replicates (if available) per site were analyzed. Only workers were taken into account to eliminate the effect of inter-caste differences in isotopic values, which could bias cross-species comparisons [35,36]. In addition, workers are the caste that does the foraging and feeding within colonies. Five replicates of soil samples were collected from the top horizon (0-15 cm) at each site-one at each study plot-resulting in 15 samples in total for all study sites. Soil samples were cooled and directly dried after each field trip, and then sealed in plastic bags. Additionally, litter samples were collected on the ground; three samples were taken per study plot, including one at the start, in the middle, and at the end of each study plot, resulting in 15 replicates per site (45 in total). They included leaves, twigs, and dead wood. Like the soil samples, they were dried and kept cool, prior to analysis. Soil samples were collected in each plot following the protocol by Pansu and Gautheyrou [47] and Osorio [48]. At a depth of 15 cm and a distance of 1 m parallel to each belt transect, three replicate soil samples were taken along the transect belt (one at the start, in the middle, and at the end of a belt transect), resulting in a total of 45 samples (3 sites × 5 belt transect × 3 replicates). Samples were prepared according to a protocol of the Centre for Stable Isotope Research and Analysis (KOSI) at the University of Göttingen (Germany). In short, all samples were dried at 60 • C for 24 h. Stones and gravel were removed before crushing samples and grounding them into fine powder. For soil and litter between 0.4 mg and 1.0 mg, one whole termite worker with gut was weighted, transferred into tin zinc capsules (HekaTech GmbH ® , Wegberg, Germany), and sent to KOSI.
Carbon and nitrogen stable isotope ratios were measured on an elemental analyzer (NA 1500, Fisons-Instruments, Rodano, Milan, Italy) and an isotope ratio mass spectrometer (Delta V Plus, Thermo Fisher Scientific, Bremen, Germany). Stable isotope ratios were expressed using the delta (δ) notation in according to: where R sample is the isotopic ratio of the sample ( 13 C/ 12 C or 15 N/ 14 N), R standard is the isotopic ratio of the international standard, and X is the respective element ( 13 C or 15 N). For 13 C V-PDB and 15 N, atmospheric nitrogen was used as the standard. Acetanilide (C 8 H 9 NO, Merck, Darmstadt, Germany) was used for internal calibration. For the amount of animal tissue analyzed per sample, precision of the measurement was about 0.1 for 13 C and 0.2 for 15 N. We calculated the mean and the standard deviation (SD) of all samples for each site.

Phylogenetic Community Analyses
The species pool of the three study sites comprises 32 species that have been morphologically and genetically identified (Table A1: GenBank accession numbers MH09082-MH090914 and KU510330, KX267100, KX267099, KX267098, KX267095, KX267092) [26]. As the input tree for the phylogenetic community structure analyses, we used the combined COII, 12S, and 16S nucleotide sequences and performed a Bayesian approach using MrBayes 3.2.1. [49]. We pruned the tree prior to analysis to include only species of the regional species pool and only one representative per species in the tree ( Figure A1, Table A2).
A commonly used index to quantify the phylogenetic structure of a local community is the Net Relatedness Index (NRI) [5]. It measures whether locally co-occurring species are phylogenetically more/less closely related than expected by chance. It uses phylogenetic branch length to measure the distance between each sample to every other terminal sample in the phylogenetic tree, hence the degree of overall clustering. It is calculated as the difference between the mean phylogenetic distance (MPD) of the tested local community (i.e., each study plots) and the MPD of the regional community (i.e., all 32 species identified for dry forests in this region), divided by the standard deviation of the latter. NRI values close to zero indicate random community assembly, which may imply that neutral processes are important in structuring communities. Large positive values reflect phylogenetic clustering of co-occurring species (i.e., co-occurring species are more related than expected by chance), whereas low negative values point to over-dispersion (i.e., co-occurring species are less related than expected by chance) [1].
Depending on whether niche-relevant traits, such as the feeding niche, are evolutionary labile or conserved, the NRI values can hind at different assembly processes [3]: For instance, conserved traits and over-dispersion can indicate that interspecific competition plays an important role in structuring communities. We analyzed the phylogenetic community structure with PHYLOCOM 4.2 [5]. As the input tree, we used the Bayesian inference tree in combination with abundance data for all species. We conducted two analyses, including one testing the local assemblages against the regional species pool (all species found during this study) and one site specific analysis in which we tested the local assemblages against the species occurring at a specific study site. We tested whether our data significantly deviated from null models using the independent swap algorithm on occurrence data [50]. This algorithm creates swapped versions of the sample/species matrix while constraining row (species) and column (occurrence) totals to match the original matrix. We used two-tailed significant rank tests as suggested by Webb et al. [5] to determine if observed values differed significantly from the null model (e.g., with 9999 randomizations, rank values equal or higher than 9750 or equal or lower than 250 are statistically significant at p = 0.05).

Mapping Food Niche Traits on Phylogeny
In order to interpret the results of the phylogenetic community analysis, it is necessary to know whether the studied food niche traits were phylogenetically conserved or labile. To determine this, we conducted two analyses, including one with feeding groups and one for isotope signatures. We used (i) the feeding groups and (ii) the mean of the δ 13 C and δ 15 N values calculated over all collection sites for each species as character states and the phylogenetic tree from Casalla & Korb [26] (which was inferred from molecular sequence data) as the input, and performed ancestral state reconstruction (ASR). For inferring ancestral states, we used Mesquite version 3.04 [51,52], in particular, the module 'Parsimonious Ancestral States: 'Parsimony unordered' for the categorical feeding group data and 'Parsimony Squared' for quantitative isotope data.

Phylogenetic Community Structure
Overall, the NRI values which measured the phylogenetic structure of the termite assemblages across the regionals species pool ranged from −0.82 to 2.45 (Table A3a). NRI values did not differ significantly from random expectation, except for one site in Colosó (Colosó 5), which showed significant signs of phylogenetic clustering (Table A3a). Both sites, Colosó and Ceibal, had significantly higher NRI values than Tayrona, where species were more phylogenetically overdispersed ( Figure 2, Table A3a). At the study plot level, NRI values did not correlate with species richness (Pearson correlation r = 0.339, p = 0.217). However, NRI values significantly increased with rainfall at a study plot (r = 0.862, p < 0.001; Figure 3a) and its elevation (r = 0.626, p = 0.012; Figure 3b). Contrarily, NRI decreased with temperature (r = −0.648, p = 0.009; Figure 3c). occurrence data [50]. This algorithm creates swapped versions of the sample/species matrix while constraining row (species) and column (occurrence) totals to match the original matrix. We used twotailed significant rank tests as suggested by Webb et al. [5] to determine if observed values differed significantly from the null model (e.g., with 9999 randomizations, rank values equal or higher than 9750 or equal or lower than 250 are statistically significant at p = 0.05).

Mapping Food Niche Traits on Phylogeny
In order to interpret the results of the phylogenetic community analysis, it is necessary to know whether the studied food niche traits were phylogenetically conserved or labile. To determine this, we conducted two analyses, including one with feeding groups and one for isotope signatures. We used (i) the feeding groups and (ii) the mean of the δ 13 C and δ 15 N values calculated over all collection sites for each species as character states and the phylogenetic tree from Casalla & Korb [26] (which was inferred from molecular sequence data) as the input, and performed ancestral state reconstruction (ASR). For inferring ancestral states, we used Mesquite version 3.04 [51,52], in particular, the module 'Parsimonious Ancestral States: 'Parsimony unordered' for the categorical feeding group data and 'Parsimony Squared' for quantitative isotope data.

Phylogenetic Community Structure
Overall, the NRI values which measured the phylogenetic structure of the termite assemblages across the regionals species pool ranged from −0.82 to 2.45 (Table A3a). NRI values did not differ significantly from random expectation, except for one site in Colosó (Colosó 5), which showed significant signs of phylogenetic clustering (Table A3a). Both sites, Colosó and Ceibal, had significantly higher NRI values than Tayrona, where species were more phylogenetically overdispersed ( Figure 2, Table A3a). At the study plot level, NRI values did not correlate with species richness (Pearson correlation r = 0.339, p = 0.217). However, NRI values significantly increased with rainfall at a study plot (r = 0.862, p < 0.001; Figure 3a) and its elevation (r = 0.626, p = 0.012; Figure 3b). Contrarily, NRI decreased with temperature (r = −0.648, p = 0.009; Figure 3c).
A mixed effect model of NRI and the three abiotic variables (rainfall, elevation, and temperature), using site as the random factor, showed significant effects for temperature only (Table 1). Mixed effects using fine-scale were not significant (p > 0.082, Table A3b).

Isotopes Stable Analyses
δ 15 N values of termites were highly variable, ranging from −1.6 to 17.8 (mean 6.0 +/− 1 SD 3.8, Figure A2, Table A4). Both Colosó and Tayrona differed significantly from Ceibal (F 2, 12 = 4.78; p = 0.03, Figure 4a, Table A5). Additionally, the δ 13 C values of termites were highly variable (mean  Table A6). In addition, the δ 15 N and δ 13 C signatures for litter samples were also significantly lower at these two sites than at the There were also significant differences in isotope signatures between feeding groups. δ 15 N values differed significantly between all groups (F1, 212 = 23.80, p < 0.001) (Figure 5a). Species from feeding group I had significantly lower δ 15 N values than those of feeding group II (p = 0.047), III (p < 0.001), and IV (p < 0.001). Termites from feeding group II had significantly lower δ 15 N values than those of the other feeding groups, III (p < 0.001) and IV (p < 0.001). Termites from feeding group III had significantly lower from feeding group IV (p < 0.001). Thus, there was a gradual increase of δ 15 N over the feeding groups.
For the δ 13 C, there were less strong differences between groups. Only feeding group II had significantly lower values than all other feeding groups (mean: −28.2‰ +/− 1 SD 1.0; F1, 212 = 9.34, p = 0.003, Figure 5b). There were also significant differences in isotope signatures between feeding groups. δ 15 N values differed significantly between all groups (F 1, 212 = 23.80, p < 0.001) (Figure 5a). Species from feeding group I had significantly lower δ 15 N values than those of feeding group II (p = 0.047), III (p < 0.001), and IV (p < 0.001). Termites from feeding group II had significantly lower δ 15 N values than those of the other feeding groups, III (p < 0.001) and IV (p < 0.001). Termites from feeding group III had significantly lower from feeding group IV (p < 0.001). Thus, there was a gradual increase of δ 15 N over the feeding groups.

Mapping Food Niche Traits on Phylogeny
Our analyses showed that food niches, measured as feeding group membership and δ 15 N and δ 13 C signatures, are phylogenetically conserved traits in the studied species. Closely related species share the same feeding group (right part of Figure 7a,b). Among the studied termites, group IV soil feeders evolved only once from group II plant litter feeders or group I wood feeders. Interestingly, group IV soil feeders do not seem to have evolved from group III soil feeders (and vice versa). At the fine-scale of the δ 15 N and δ 13 C signatures, the δ 15 N signal reflects the feeding group pattern well, except for a few species, such as Incisitermes schwarzi and Termes sp.1 (Figure 7a, left part). Thus, the δ 15 N signature has a strong phylogenetic signal, with closely related species sharing similar signatures (Figure 7a). The δ 13 C signatures are also phylogenetically conserved, but their pattern does not reflect that of the feeding groups (Figure 7b, right part). The Kalotermitidae (feeding group I) (especially Cryptotermes) had the highest δ 13 C values, while Nasutitermitinae and Microcerotermes (both feeding group II, but independent transitions) had the lowest values. For δ 13 C, the subterranean Rhinotermitidae had the highest values (−26.4 +/− 1 SD 1.0), which were significantly higher than for species from the Termitidae (−27.1 +/− 1 SD 1.2, Figure 6b).

Mapping Food Niche Traits on Phylogeny
Our analyses showed that food niches, measured as feeding group membership and δ 15 N and δ 13 C signatures, are phylogenetically conserved traits in the studied species. Closely related species share the same feeding group (right part of Figure 7a,b). Among the studied termites, group IV soil feeders evolved only once from group II plant litter feeders or group I wood feeders. Interestingly, group IV soil feeders do not seem to have evolved from group III soil feeders (and vice versa). At the fine-scale of the δ 15 N and δ 13 C signatures, the δ 15 N signal reflects the feeding group pattern well, except for a few species, such as Incisitermes schwarzi and Termes sp.1 (Figure 7a, left part). Thus, the δ 15 N signature has a strong phylogenetic signal, with closely related species sharing similar signatures ( Figure 7a). The δ 13 C signatures are also phylogenetically conserved, but their pattern does not reflect that of the feeding groups (Figure 7b, right part). The Kalotermitidae (feeding group I) (especially Cryptotermes) had the highest δ 13 C values, while Nasutitermitinae and Microcerotermes (both feeding group II, but independent transitions) had the lowest values.   [25]; Dark blue represents feeding group I, light blue represents feeding group II (down-side tree), green represents feeding group III, and red represents feeding group IV. Ancestral states represented by colors at the nodes of the phylogeny observed at the tips are circled at each node. Families: K: Kalotermitidae, R: Rhinotermitidae, and T: Termitidae. Subfamilies: Ap: Apicotermitinae, Na: Nasutitermitinae, Te: Termitinae, and Sy: Syntermitinae.  [25]; Dark blue represents feeding group I, light blue represents feeding group II (down-side tree), green represents feeding group III, and red represents feeding group IV. Ancestral states represented by colors at the nodes of the phylogeny observed at the tips are circled at each node. Families: K: Kalotermitidae, R: Rhinotermitidae, and T: Termitidae. Subfamilies: Ap: Apicotermitinae, Na: Nasutitermitinae, Te: Termitinae, and Sy: Syntermitinae.

Discussion
Our results imply that mechanisms which structure the termite assemblages differ between sites, with interspecific competition being more important at drier and warmer, lower-altitude plots (Figures 2 and 3). This is in line with a hypothesis of food-limitation becoming important in such areas.

Mechanisms Structuring Termite Assemblages
Inferred from the phylogenetic community analyses, the assembly processes in the studied Colombian dry forests seem to differ between sites. The driest and lowest elevation site, Tayrona, had termite assemblages that were phylogenetically more overdispersed than those of the other sites ( Figure 2). Overall, phylogenetic overdispersion correlated negatively with rainfall and elevation of study plots and positively with temperature (or vice versa, phylogenetic clustering increased with rainfall and elevation, but decreased with temperature) (Figure 3). The mixed model analyses including all three environmental analyses together, revealed that rainfall had the strongest and only significant effect (Table 1). That the results are non-significant when using the NRI values generated with the site-specific species pools as the reference supports the conclusion that the pattern is mainly driven by abiotic differences between study sites and not by differences between study plots within a site. To interpret this pattern ecologically requires knowledge of whether niche traits are phylogenetically conserved or labile. For conserved traits (i.e., closely related species share traits), overdispersion implies that species which share the same niche traits are less likely to co-exist than species that differ in these traits. Our phylogenetic trait mapping analyses showed that the food niche traits are phylogenetically conserved for the studied termite assemblages. This indicates that interspecific competition is more important in Tayrona, and in general, at the drier and warmer low-elevation study plots, than at the more humid and slightly colder high-elevation plots. In line with this, there are fewer termite species at the Tayrona plots [26]. In tropical dry forests, rainfall rather than temperature limits vegetation growth [53], with the latter still being optimal for plant growth at all sites (mean 26.5-27.7 • C). In line, the mixed model analyses only showed a significant effect for rainfall. Thus, there is lusher vegetation with higher biomass production in Colosó [54] and Ceibal [55,56] than in Tayrona [41], where Euphorbiaceae and Cactaceae are common. This supports the hypothesis that food is a limiting resource over which termites compete at Tayrona. Furthermore, other studies found evidence that food can be a limiting resource (i.e., dead plant material) for termites [57], which can lead to intra-and interspecific competition [57][58][59][60][61]. In a West African savannah, annual fires reduce the availability of dead plant material, so that the addition of dried grass after the fires leads to an increase in the number of sexual produced by the dominant mound building termite Macrotermes bellicosus [57].
This does not seem to be the case at the two other sites. NRI values close to zero imply that assemblages at Colosó and Ceibal do not differ much from random associations. The fact that the NRI values at the plot scale did not differ significantly from random expectation might be due to the low numbers of species per plot. Note, that the species number is the sample size for these tests and that small sample sizes are associated with low statistical power, hence making it unlikely to detect significant effects. One plot in Colosó even showed signs of phylogenetic clustering, implying that environmental variables select for a certain subset of termites at this specific plot, which was characterized by huge trees with a dense canopy and high humidity. Four Nasutitermes species and four species of the Anoplotermes-group co-existed in this plot. Several studies from Neotropical rain forests have also shown that these closely related species commonly co-occur [62][63][64].
How wide-spread are the implicated structuring mechanisms in termites? Comparable studies are rare. Most research on termite communities has concentrated on describing local or regional assemblages and testing associations between termite diversity and variables such as fire, disturbance, or elevation gradients [65][66][67][68]. The few studies that have addressed community processes in more detail in termites have often concentrated on a subset of species. They found evidence for interspecific competition in structuring assemblies at the local scale [32,44,58,69,70]. The only studies directly comparable to our study come from Africa. There phylogenetic community analyses imply random processes [49], but also evidence for interspecific competition [12,71] and environmental filtering [13], depending on the study site, disturbance regime, and presence of a dominant mound building termite species. Thus, we currently cannot derive any general conclusions and more similar studies are needed.

Food Niche: Isotopes and Termite Feeding Groups
In our study, also the isotope analyses, which admittedly included the termites' gut, of the litter support the notion of vegetation differences in Tayrona as the δ 15 N and δ 13 C values were significantly higher at this site. Interestingly, the soil signatures did not differ between sites. For δ 13 C, but not δ 15 N, the shift in the litter signatures between sites is reflected in the termites' isotope signal ( Figure 5). Differences in δ 13 C mainly reflect varying proportions of C 3 and C 4 plants in the food of termites, while differences in δ 15 N indicate the diverse proportions of variably humified food resources [32,34,36]. The higher δ 13 C values at Tayrona reflect the presence of more grasses and especially abundant Euphorbiaceae [41], which are all C 4 plants. In addition, due to the proximity of the sea and the associated salinity-and water-stress, also C 3 plants have higher δ 13 C values [72,73].
In general, our isotope signatures were similar to those found for termites in other forests [31,32]. However, a study in an African savanna with many fungus-growing termites revealed higher δ 13 C values at the upper range and a lower δ 15 N signature at the lower range, reflecting a higher proportion of C4 grasses in the habitat and a broader food niche spectrum of fungus-growing termites [33,37].
Reflecting a humification gradient, the four commonly recognized feeding groups identified by Donovan et al. [27] can generally be distinguished by isotope signatures, especially δ 15 N [34,37]. Nevertheless, there can be limitations, as the rain forest and savanna revealed, which did not recover a discontinuity in δ 15 N values between group I and II or between III humus-and group IV soil-feeders [74]. Our current study separated all feeding groups for δ 15 N signatures (Figures 5a, 6a and 7a), supporting Donovan´s feeding group concept (Figure 7b). However, for δ 13 C, no such gradient was revealed and only feeding type II had δ 13 C values that were lower and differed significantly from all other feeding types (Figures 5b, 6b and 7b). This implies that isotope studies are required to reliably determine the food niches of termites.

Rainfall
Elevation Temperature Estimates of the mixed effect model between NRI (using the site-specific termite pools) and the three abiotic variables rainfall, elevation, and temperature. The mixed model was insignificant (p > 0.087).  Figure A2. Scatterplot between δ 15 N and δ 13 C. Symbols: Orange square: feeding group I; green circle: feeding group II; purple circle: feeding group III; black circle: feeding group IV.  . Scatterplot between δ 15 N and δ 13 C. Symbols: Orange square: feeding group I; green circle: feeding group II; purple circle: feeding group III; black circle: feeding group IV.

Variable
Insects 2019, 10 18 of 26 Figure A2. Scatterplot between δ 15 N and δ 13 C. Symbols: Orange square: feeding group I; green circle: feeding group II; purple circle: feeding group III; black circle: feeding group IV.