Background: Recent increases in understanding the ecological and evolutionary roles of microbial communities has underscored their importance for their hosts’ biology. Yet, little is known about gut microbiota dynamics during early stages of ecological diversification and speciation. We sequenced the V4 region of the 16s rRNA gene to study the gut microbiota of extremely young adaptive radiations of Nicaraguan Midas cichlid fish (Amphilophus cf. citrinellus) from two crater lakes to test the hypothesis that parallel divergence in trophic ecology is associated with parallel changes of the gut microbiota.

Results: Bacterial communities of the water and guts were highly distinct, indicating that the gut microbiota is shaped by host-specific factors. Across individuals of the same crater lake, differentiation in trophic ecology was 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 two crater lakes, we could not detect evidence for parallel changes of the gut microbiota associated with trophic ecology.

Conclusions: Similar cases of non-parallelism have been observed in other recently diverged fish species and might be explained by a lack of clearly differentiated niches during early stages of ecological diversification.

Specimens of the Amphilophus cf. citrinellus species complex were caught during field trips to Nicaragua in 2014 and 2015. 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 water samples were collected along the shores of the four lakes in 2018. Briefly, for each replicate 500 ml of lake water were filtered through a cellulose nitrate filter (ø 25 mm, pore size 1 µm) and filters were stored in Longmire’s solution  at -20°C until DNA extraction. DNA from midgut tissue was extracted using the commercial QIAamp DNA Stool Mini Kit according to the manufacturer’s protocol (Qiagen, Hilden, Germany). DNA from water samples was extracted from cellulose nitrate filters using a QIAGEN DNeasy Blood & Tissue kit. All DNA extractions and PCR amplifications were performed under sterile conditions in a laminar flow hood to minimize contamination risk. Fish gut and water samples were separately pooled in an equimolar manner and size selection was performed on a Pippin Prep device (Sage Science, Beverly, MA). The quality of the pooled libraries was assessed using a Bioanalyzer 2100 (Agilent Technologies, Waldbronn, Germany). Both libraries were paired-end sequenced, each in one lane of the Illumina flow cell. For the fish guts, we sequenced 2x250 bp on an Illumina HiSeq 2500 platform at TUCF Genomics (Tufts University). Water samples were sequenced 2x150 bp on an Illumina HiSeq X-ten at BGI Genomics.

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 water samples, respectively. Illumina adapters were removed and reads were trimmed with Trimmomatic v0.36. As there was no overlap between forward and reverse reads for water samples and the sequence quality of forward reads was higher, we used 135 bp of the forward reads for all analyses. The demultiplexed and trimmed reads were imported into the open-source bioinformatics pipeline Quantitative Insights Into Microbial Ecology [QIIME2] to analyze microbial communities of fish guts and water samples. Briefly, sequence quality control was done with the QIIME2 plugin deblur. A phylogenetic tree of bacterial taxa was produced with FastTree 2.1.3. Bacterial species richness and bacterial community composition [weighted UniFrac, unweighted UniFrac, Bray-Curtis dissimilarity] were calculated and taxonomy was assigned using vsearch against the SILVA 132 ribosomal RNA (rRNA) databases at a 97% similarity threshold. Statistical analyses were performed in R v3.2.3.