Parallel and non-parallel changes of the gut microbiota during trophic diversification in repeated young adaptive radiations of sympatric cichlid fish

Background Recent increases in understanding the ecological and evolutionary roles of microbial communities have underscored the importance of their hosts’ biology. Yet, little is known about gut microbiota dynamics during the early stages of ecological diversification and speciation. We sequenced the V4 region of the 16s rRNA gene to study the gut microbiota of Nicaraguan Midas cichlid fish (Amphilophus cf. citrinellus). Specifically, we tested the hypothesis that parallel divergence in trophic ecology in extremely young adaptive radiations from two crater lakes is associated with parallel changes of their gut microbiota. Results Bacterial communities of fish guts and lake water were highly distinct, indicating that the gut microbiota is shaped by host-specific factors. Among individuals of the same crater lake, differentiation in trophic ecology was weakly associated with gut microbiota differentiation, suggesting that diet, to some extent, affects the gut microbiota. However, differences in trophic ecology were much more pronounced across than within species whereas similar patterns were not observed for taxonomic and functional differences of the gut microbiota. Across the two crater lakes, we could not detect conclusive evidence for parallel changes of the gut microbiota associated with trophic ecology. Conclusions A lack of clearly differentiated niches during the early stages of ecological diversification might result in non-parallel changes of gut microbial communities, as observed in our study system as well as in other recently diverged fish species. Video Abstract


Introduction
The importance of microorganisms for many aspects of their hosts' biology is increasingly recognized for a wide range of animals, from insects to mammals (1)(2)(3). The gut microbiota is a complex and dynamic community that is fundamental for physiological processes, such as regulation of the immune system (4) and nutrient metabolism (5). Further, the significance of microbes in animal evolution has become increasingly appreciated (1,6,7). In some cases, divergence of the gut microbiota appears to be strongly correlated with their host's phylogeny and genetic divergence (8,9). These findings are supported by the fact that host genetics, together with environmental effects such as diet, contributes to shaping and maintaining gut microbiota composition (10)(11)(12). However, open questions remain on how closely the gut microbiota matches the biology of its host during ecological diversification and speciation or whether the composition of the gut microbiota could even be predicted based on the ecology of its host. These important questions can best be addressed in a setting where evolution repeated itself, i.e., in pairs of species that evolved in parallel. Cases of parallel evolution of host species associated with divergence in trophic ecology and habitat use allow us to ask whether the gut microbiota also changes in a predictable and parallel manner. This question has been addressed in several fish species covering a wide range of divergence times but results have been inconsistent. African cichlids from two old adaptive radiations of Barombi Mbo (0.5 -1 myr) and Tanganyika (9 -12 myr) show parallel changes of the gut microbiota associated with host diet (13). Yet, studies on lineages that diverged more recently like whitefish and Trinidadian guppies did not find evidence for parallelism (14,15).
Here, we asked whether evolutionary divergence in host species' trophic ecology predicts the composition of the gut microbiota of Nicaraguan Midas cichlids (Amphilophus cf. citrinellus), which represent very young adaptive radiations that, notably, occurred in parallel and sympatrically. Currently, there are 13 described species of Midas cichlids (16)(17)(18)(19) and their distribution results from independent colonization events from two older great lakes (Lakes Managua and Nicaragua) that are approximately 500,000 years old (17,20) into several crater lakes (calderas of inactive volcanoes; all crater lakes are between 1,000 and 23,000 years old) (21). Crater lake Midas cichlids differ from their source populations of the great lakes in traits such as body shape and visual sensitivity (22)(23)(24)(25). The colonization events of crater lakes Apoyo (colonized from Lake Nicaragua) and Xiloá (colonized from Lake Managua) are estimated to have occurred as recently as 1,700 and 1,300 generations ago, respectively (18,26). Within these two crater lakes, multiple species of Midas cichlids evolved in sympatry during these extremely short time spans, hence, ten species are endemic to Apoyo (six) and Xiloá (four) (18,27). Notably, one slenderbodied, limnetic species (A. zaliosus in Apoyo and A. sagittae in Xiloá) independently evolved in each of the two crater lakes. These elongated limnetic species are not found in the great lakes and inhabit the open water zone that is exclusive to the deep crater lakes. Limnetic species differ distinctly in body shape from several deep-bodied benthic species in their respective crater lakes (18,22,25) and feed at a higher trophic level (22) (Fig. 1). Previously, we have shown that gut microbiotas differ between a benthiclimnetic species pair from Apoyo (28).
The extraordinary system of crater lake Midas cichlids is a promising model to elucidate to what extent trophic ecology is mirrored by repeated and parallel changes of the gut microbiota across two very young adaptive radiations. First, we investigated trophic position (stable isotope ratios of carbon and nitrogen) and the gut microbiota of Midas cichlids from the two source lakes, the great lakes Managua and Nicaragua, four species from crater lake Apoyo and three species from crater lake Xiloá. in particular, we tested the hypotheses that (i) species from distinct lakes differ in their gut microbiota but also from the bacterial communities of their natural environments and (ii) repeated adaptation to different trophic niches is associated with parallel changes of the gut microbiota across the two crater lake radiations.

Sample collection
Specimens of the Amphilophus cf. citrinellus species complex were caught during field trips to Nicaragua . We collected A. citrinellus populations from the two great lakes, Lake Nicaragua (n = 10) and Lake Managua (n = 10).
From two crater lakes, we collected elongated limnetic species, A. sagittae from Xiloá (n = 20) and A. zaliosus from Apoyo (n = 19), as well as several deep-bodied benthic species, A. amarillo (n = 17) and A. xilaoensis (n = 16) from Xiloá, A. astorquii (n = 19), A. chancho (n = 11) and A. globosus (n = 5) from Apoyo ( Fig. 1). All specimens were sacrificed by applying an overdose of MS-222. Then, whole guts were dissected, cleaned and stored in absolute EtOH at -20°C until DNA extraction. Muscle tissue of the same specimens was collected in absolute EtOH and stored at -20°C for stable isotope analyses. Four technical replicates of environmental DNA (eDNA) were collected along the shores of the four lakes in 2018. Briefly, 500 ml of lake water were filtered through a cellulose nitrate filter (ø 47 mm, pore size 1 µm), stored in Longmire's solution (29) at -20°C until DNA extraction.

Trophic analysis
Stable isotope ratios of carbon (δ 13 C) and nitrogen (δ 15 N) were determined based on muscle tissue of the same fish used for gut microbiota analyses. Dried and powdered samples (0.6 mg) were loaded into tin capsules and combusted in a vario Micro cube elemental analyzer (Elementar Analysensysteme, Germany). The resulting gases were fed via gas chromatography into the inlet of a Micromass Isoprime Isotope Ratio Mass Spectrometer (Isoprime, Cheadle Hulme, UK). Two sulfanilamides (Iso-prime internal standards) and two Casein standards were used. Internal laboratory standards indicated measurement errors (SD) of ± 0.03‰ for δ 13 C and 0.12‰ for δ 15 N. Isotopic values are reported in δ-notation in parts per thousand deviations (‰) relative to international standards for carbon (Pee Dee Belemnite, PDB) and nitrogen (atmospheric N2, AIR) according to the following equation: (‰) = 1000 * ( − 1)

Library preparation & Illumina Sequencing
DNA from fish guts was extracted from approximately 50-100 mg of tissue from the medial gut section

Gut microbiota analysis
We obtained a total of 62,728,287 (median: 238,073 reads/specimen) and 111,949,556 (median: 6,825,739) raw sequencing reads that could be unambiguously assigned to a specific sample for fish guts and eDNA, respectively. Illumina adapters were removed and reads were trimmed with Trimmomatic v0.36 (30). As there was no overlap between forward and reverse reads for eDNA samples and the sequence quality of forward reads was higher, we used 135 bp of the forward reads for all analyses. The databases at a 97% similarity threshold (36). The weighted UniFrac distance matrix was visualized with principal coordinate analyses. Since most diversity metrics are sensitive to differences in numbers of reads per sample, we chose 20,000 sequencing reads, the approximate number of reads for the sample with the lowest sequencing depth, as our sampling depth for all further analyses. To determine whether this sampling depth was appropriate to capture a large proportion of the microbial diversity for each sample (bacterial species richness measured as observed OTUs), we rarefied our data (Fig. S1). This analysis confirmed that a large proportion of the gut microbial diversity is already captured at a sequencing depth of 20,000 reads/sample (Fig. S1). Non-parametric Wilcoxon rank-sum tests were used for pairwise comparisons (37) and Kruskal-Wallis tests for comparisons across multiple groups, as implemented in the R stats package (38

Diet differentiation among Midas cichlids
In order to obtain information on trophic ecology of all studied Midas cichlid species, we measured stable isotope ratios of carbon (δ 13 C) and nitrogen (δ 15 N). Overall, Midas cichlids from different environments significantly differed in δ 13 C (Kruskal-Wallis test, P < 0.001) and δ 15 N (P < 0.001; Fig. 1). Within each of the crater lakes, stable isotope ratios differed among species (δ 13 C: PApoyo < 0.001, PXiloá = 0.001; δ 15 N: P < 0.001 for both lakes). These results illustrate that sympatric species of Midas cichlids feed on different carbon sources and further suggests that they also occupy distinctive trophic levels based on nitrogen values (22).
As predicted based on diet and inferred trophic niche (22,27), the limnetic species had the highest nitrogen value in both crater lakes (Fig. 1). Benthic species occupied distinct trophic niches that were generally at lower trophic levels (A. globosus in Apoyo, A. amarillo in Xiloá) than the limnetic species. Yet, one benthic species largely overlapped with the limnetic species in each crater lake (A. chancho with the limnetic A. zaliosus in Apoyo, A. xiloaensis with the limnetic A. sagittae in Xiloá). The benthic A. astorquii from Apoyo was highly variable in carbon and nitrogen signatures and largely overlapped with the other species (Fig. 1).

Gut microbiota differentiation across lakes
Bacterial community composition significantly differed between lake water and fish guts (weighted UniFrac; adonis, P = 0.001; Fig. 2), emphasizing that the gut microbiota not merely represents the microbial community of the natural environment. In the water samples, Cyanobacteria (9.9 -21.9%), Planctomycetes (10 -23%) and Actinobacteria (9 -25.1%) constituted a large proportion of microbial communities whereas these groups where much less abundant in the gut microbiota of Midas cichlids (Fig. 3A). In contrast, the OTUs) (Wilcoxon rank-sum test, P < 0.001; Fig. 3A) and among the water samples, great lakes showed a higher bacterial species richness than crater lakes (P < 0.001). Among Midas cichlids, bacterial community composition differed not only across lakes (P = 0.001) but also between environment type (great lakes vs. crater lakes; P = 0.007). Bacterial species richness was lower in great lake Midas cichlids (P < 0.001), however, this pattern disappeared when A. citrinellus from Lake Managua was removed from the analysis (P = 0.5124). These results clearly show that bacterial species richness is largely constant across Midas cichlids from different environments, except for A. citrinellus from Lake Managua that showed strongly reduced bacterial diversity. Next, we investigated whether gut microbiota divergence within crater lakes is associated with the observed differences in stable isotope ratios.

Association between trophic ecology and gut microbiota in crater lake adaptive radiations
The two crater lakes are each inhabited by several endemic Midas cichlid species that substantially differ in their morphology, ecological niche, and diet (Fig. 1). Yet, there were no significant differences in bacterial species richness (observed OTUs) among sympatric species within each of the crater lakes (Wilcoxon rank-sum tests, P > 0.05 for all pairwise comparisons). Hence, we tested whether bacterial community composition (weighted UniFrac) varied among species of the two parallel crater lake radiations of Apoyo and Xiloá (Fig. S2).
There was overall differentiation in the taxonomic composition of bacterial communities among sympatric species in both crater lakes (PApoyo = 0.001, PXiloá = 0.007). From the taxonomic composition of the gut microbiota, one can infer the functional bacterial metagenome by predicting the abundance of genes involved in metabolic pathways (40). Predicted functional bacterial metagenomes significantly differed across species only in crater lake Xiloá (P = 0.006). Since there is pronounced variation in trophic ecology of both adaptive radiations (Fig. 1), we calculated pairwise distances of trophic ecology (δ 15 N and δ 13 C scores) and correlated these with pairwise distances in bacterial community composition among all individuals within each crater lake (Fig. 4). For carbon, we found a significant positive correlation with bacterial community composition in Apoyo (Pearson's product-moment correlation, r = 0.092, P < 0.001; Distinguishing between intra-and interspecific comparisons of pairwise distances revealed that differentiation in diet is highly distinct among than within species for carbon (Wilcoxon rank-sum test, P < 0.001 for both lakes; Fig. 5A) and nitrogen (P < 0.001 for both lakes; Fig. 5B). In contrast, taxonomic differentiation of the gut microbiota (PApoyo = 0.069, PXiloá = 0.004; Fig. 5C) and the predicted functional metagenome (PApoyo = 0.824, PXiloá = 0.047; Fig. 5D) showed more similar (albeit significantly different in Xiloá) levels between intra-and interspecific comparisons. This clearly illustrates that pronounced differentiation in diet across species is not reflected by equivalent changes of the gut microbiota within crater lakes.

Parallelism and non-parallelism of the gut microbiota in crater lake Midas cichlids
Next, we investigated whether the repeated evolution of sympatric species differing in trophic ecology in the crater lakes is associated with parallel changes of the gut microbiota. As δ 15 N values were consistently higher in crater lake Xiloá compared to crater lake Apoyo, we performed z-score normalization of the data to allow comparisons across crater lakes by inferring trophic position (normalized δ 15 N) and littoral carbon usage (normalized δ 13 C) ( Fig. S3; 42). these results indicate that overall gut microbiota differentiation did not occur in parallel with divergence in Midas cichlids' trophic ecology across the two crater lakes.

Discussion
Numerous studies on diverse vertebrate species have convincingly demonstrated that the composition of bacterial gut communities is affected by diet (13,(43)(44)(45). What remains largely unknown are gut microbiota dynamics during the host's adaptation to novel food sources, particular during early stages of species divergence. To address this question, we studied trophic ecology and the gut microbiota of repeated Nicaraguan Midas cichlid crater lake radiations, a model system for rapid ecological diversification and speciation (17,18,27). We asked whether the parallel evolution of trophic diversification is associated with respective changes of the gut microbiota among sympatric crater lake species (i.e., are the gut microbiotas of ecologically similar species that independently evolved in two crater lakes more similar to each other than they are to their closest relatives of the same lake?). Our results clearly demonstrate that among individuals of the same crater lake differentiation in trophic ecology and the gut microbiota is associated, emphasizing the importance of diet in shaping the gut microbiota. However, we did not find evidence for parallel changes of the gut microbiota across crater lakes, suggesting that diet affected Midas cichlids' gut microbiota differently in these lakes. Moreover, we found that interspecific variation in trophic ecology (measured as stable isotope ratios of carbon and nitrogen) is significantly higher than intraspecific variation, a pattern that was much less pronounced for the gut microbiota (Fig. 4).

Gut microbiota differentiation across lakes
Comparing bacterial communities from environmental DNA (eDNA) with those harbored in fish guts clearly revealed that, although bacterial taxa are shared to some extent, the gut microbiota not merely represents the bacterial community of the natural environment but is rather controlled by the host, as has been found for other fishes (14,15,46). The dominant bacterial phyla of Midas cichlids' gut microbiota (Proteobacteria, Firmicutes, Fusobacteria and Bacteroidetes) are also found in many other freshwater fishes (14,15,47). Albeit the bacterial species richness of environmental samples strongly differed across lakes (Fig. 3B), the gut microbiota of Midas cichlids, except for A. citrinellus from Lake Managua, showed constant levels for this measure. This indicates that the diversity of bacterial species in the gut might be constrained by the host and stabilized at a given level, as predicted by the holobiont theory (48). In A.
citrinellus from Lake Managua, bacterial species richness was by far the lowest among all populations (Fig.   3B) and was also lower in eDNA from Lake Managua compared to Lake Nicaragua. The city of Managua, Nicaragua's capital with a population of more than 2 million, is located on the shore of Lake Managua and for decades, domestic and industrial waste water has been disposed into the lake (49). As a result, concentrations of mercury and other toxic substances are extremely high in the lake and are also enriched in fishes (49,50). Mercury levels have been shown to be correlated with δ 15 N values (51,52) and Midas cichlids from Lake Managua showed the highest δ 15 N values (Fig. 1), in agreement with the observation that mercury accumulates in these fishes. Albeit speculative at this point, high levels of contamination might have decreased the bacterial species richness of Lake Managua as well as the gut microbiota diversity of Midas cichlids inhabiting this lake.

Association between trophic ecology and gut microbiota in crater lake adaptive radiations
Diet is widely known to influence the composition of the gut microbiota (13,44,45,56) but most of the previous research mainly focused on model organisms or on distantly related species with long divergence times. However, the dynamics of gut microbiota changes during early stages of ecological diversification and speciation remain largely unknown, but see (14,15). Midas cichlids represent an excellent model to study such dynamics as species from crater lakes Apoyo and Xiloá repeatedly diverged only very recently and show strong differentiation in trophic ecology (Fig. 1). Therefore, we tested whether gut microbiota differentiation is associated with trophic ecology of crater lake Midas cichlids.
There was a positive correlation between interindividual differentiation of the gut microbiota with carbon isotope signatures in both crater lakes and with nitrogen isotope signatures in crater lake Xiloá (Fig.   4). Accordingly, divergence of trophic ecology and the gut microbiota are to a certain degree associated, suggesting that adaptation to different food sources necessitated changes of the gut microbiota. However, it should be noted that a substantial amount of gut microbiota variation is not explained by trophic divergence. Differences in stable isotope ratios, reflecting the host's trophic ecology, where considerably higher across species compared to within species (Fig. 5A & B). In contrast, such differences were much less pronounced in crater lake Xiloá and absent in Apoyo for the taxonomic composition of the gut microbiota and the predicted functional bacterial metagenome ( Fig. 5C & D). This could mean that occupation of novel trophic niches might be achieved without drastically changing the overall composition of the gut microbiota, both taxonomically and functionally. Rather, subtle changes in some functionally important bacterial taxa might suffice to exploit new food sources and to allow ecological and evolutionary diversification of the hosts. Alternatively, the very recent divergence of crater lake Midas cichlids might impede clear differentiation of the gut microbiota, as discussed in more detail in the following paragraph.

Parallelism and non-parallelism of the gut microbiota in fishes
Parallel changes of the gut microbiota associated with differentiation in trophic ecology have been reported for older fish species (13), whereas other studies found no evidence for parallelism among more recently diverged populations (14,15). In the very recent Midas cichlid adaptive radiations from Apoyo and Xiloá that diverged less than 1,700 and 1,300 generations ago, respectively (18), parallel changes in diet led us to expect that a similar pattern could be found in the gut microbiota as well. However, we did not detect evidence for parallel changes of the gut microbiota, except for a significant association of the predicted functional metagenome with trophic position (measured as normalized δ 15 N). Taken together, studies across multiple groups of fishes suggest that parallel changes of the gut microbiota might only be expected on longer time scales (13). These observations can be explained by the fact that during early stages of divergence, species might occupy novel niches but diet, to some extent, still overlaps across young species or ecotypes.
In crater lake Midas cichlids, stable isotope analyses clearly showed that species largely occupy distinct niches (Fig. 1). However, stomach content analyses showed that these species mainly feed on similar food items but their relative proportions differ across species (22). Young crater lake species might be in the process of adapting to specialized ecological niches but currently they are still opportunistic generalists with a varying diet. Thus, short term changes of an individual's diet are not reflected in the stable isotope signature of muscle tissue as this represents an average of this individual's diet over a period of approximately three months (57,58). In contrast, the composition of the gut microbiota is highly variable and changes rapidly with diet (59,60). Hence, the gut microbiota rather represents a snapshot of an individual's most recently acquired food items, generating high levels of intraspecific variation. This could explain why intra-and interspecific variation of the gut microbiota is much more equal compared to stable isotope data (Fig. 5). Accordingly, high levels of dietary intraspecific variation might mask interspecific differences in trophic ecology, thereby blurring any signal of gut microbiota parallelism in recently diverged ecotypes or species. This is what we can also see in other fishes like whitefish and guppies, where the main change between ecotypes is in the relative proportion of food items (14,15).
Only after species sufficiently diverged to become trophic specialists that do not overlap in food items (e.g., carnivorous and herbivorous cichlids from Lake Tanganyika), one would expect persistent and parallel patterns of gut microbiota divergence, as seen in African cichlids (13). Thus, taking into account ecological characteristics (e.g., extent of trophic divergence) as well as the evolutionary history of host species will aid in predicting when to expect parallel changes of the gut microbiota. This, in turn, will improve our understanding of gut microbiota dynamics during early stages of ecological diversification.
To conclude, in this study we illustrate that (i) the gut microbiota of Midas cichlids from distinct environments is clearly differentiated, (ii) within crater lakes, divergence of the gut microbiota is associated with trophic ecology but such changes did not occur in parallel across crater lakes and (iii) differentiation in trophic ecology is not reflected by large rearrangements of gut bacterial communities in these recently diverged species. in Apoyo) occur, which evolved rapidly within less than 2,000 generations. White inset: Stable isotope analysis of carbon (δ 13 C) and nitrogen (δ 15 N) based on muscle tissue. δ 13 C and δ 15 N values vary significantly not only among environments, but also within each of the crater lakes where substantial variation can be observed among species.

Figures
Within each of the crater lake, variation in δ 15 N indicates differences in trophic ecology among species.    For the gut microbiota, both taxonomically and functionally, distances are much more similar and significant differences could only be seen in crater lake Xiloá (Wilcoxon rank-sum test, *P < 0.05, **P < 0.01, ***P < 0.001).