Diving into divergence: differentiation in swimming performances, physiology and gene expression between locally-adapted sympatric cichlid fishes

Sympatric speciation occurs without geographical barriers and is thought to often be driven by ecological specialization of individuals that eventually diverge genetically and phenotypically. Distinct morphologies between sympatric populations occupying different niches have been interpreted as such differentiating adaptive phenotypes, yet differences in performance and thus likely adaptiveness between them were rarely tested. Here, we investigated if divergent body shapes of two sympatric crater lake cichlid species from Nicaragua, one being a shore-associated (benthic) species while the other prefers the open water zones (limnetic), affect cruising (U crit ) and sprinting (U sprint ) swimming abilities - performances particularly relevant to their respective lifestyles. Furthermore, we investigated species differences in oxygen consumption (MO 2 ) across different swimming speeds and compare gene expression in gills and white muscle at rest and during exercise. We found a superior cruising ability in the limnetic Amphilophus zaliosus compared to the benthic A. astorquii , while sprinting was not different, suggesting that their distinct morphologies affect swimming performance. Increased cruising swimming ability in A. zaliosus was linked to a higher oxygen demand during activity (but not rest), indicating different metabolic rates during exercise - a hypothesis supported by coinciding gene expression patterns of gill transcriptomes. We identified differentially expressed genes linked to swimming physiology, regulation of swimming behaviour and oxygen intake. A combination of physiological and morphological differences may thus underlie adaptations to these species’ distinct niches. This complex ecological specialization likely resulted in morphological and physiological trade-offs that contributed to the rapid establishment and maintenance of divergence with gene flow. with gene flow along the benthic-limnetic We tested these two species for differences in swimming performances that are thought to be adaptive in their respective divergent ecological niches (benthic vs. limnetic; Supporting References). Our results indicate that these species differ in cruising swimming performance, oxygen demand and gene expression during activity. We provide empirical evidence that external body morphology together with physiological mechanisms indeed constitute ecologically relevant traits that likely represent adaptations to the limnetic and benthic habitats. These traits facilitate niche segregation and the inferred performance trade-off likely contributes to adaptive divergence further supporting the hypothesis of ecological speciation. ) differentially

Some of the best studied examples of adaptive divergence linked to ecological specialization come from cichlid fishes, a group of freshwater fish renowned for their remarkable diversity, rapid adaptive radiations and explosive rates of speciation (e.g., Henning & Meyer, 2014;Salzburger, 2009;Seehausen, 2006). Among cichlid fishes, the Nicaraguan Midas cichlid species complex (Amphilophus spp.) with its small-scale crater lakes radiations offers a suitable model to explore population divergence due to adaptation to different ecological niches in small, contained geographic settings without geographic barriers (e.g., Barluenga et al., 2006;Coyne, 2007;Kautt et al., 2018;Meyer, 1990): Midas cichlids from the much larger, but shallow and turbid Lake Nicaragua and Lake Managua repeatedly colonized a nearby chain of small, deep and clearer crater lakes within the last few thousand years. In several of these crater lakes, the colonizers adaptively diversified extremely rapidly in sympatry, resulting in several endemic, ecologically specialized species (e.g., Barluenga et al., 2006;Kautt et al., 2018;Meyer, 1990). In Lake Apoyo, the oldest of the Nicaraguan crater lakes (approximately 24,000 years, Kutterolf et al., 2007), the endemic species A. astorquii (one of the five benthic species in this lake) and A. zaliosus (the single limnetic species) emerged within a surprisingly short amount of time (Fig. 1AB, Barlow & Munsey, 1976;Barluenga et al., 2006;Kautt et al., 2016a;Stauffer et al., 2008). The body morphologies of these fishes are thought to have adaptively differentiated from the ancestral rather benthic form into a deep-bodied benthic (A. astorquii) and an elongated limnetic shape (A. zaliosus;e.g., Elmer et al., 2010;Meyer, 1990;Recknagel et al., 2014). A. astorquii and A. zaliosus also differ in their trophic ecology: stable isotopes, scuba observations, and gut and stomach content analyses Accepted Article revealed that A. zaliosus mainly feed in the limnetic habitat on prey of a higher trophic level (insects, larvae, fish) compared to A. astorquii, which is a benthic omnivore that relies on a higher proportion of plants, algae and molluscs (Fig. 1, Barluenga et al., 2006;Elmer et al., 2014;Franchini et al., 2014a). These two species also host divergent gut bacterial communities (Franchini et al., 2014a) and mate assortatively (Barlow & Munsey, 1976;Baylis, 1976;Stauffer et al., 2008). These features make Neotropical Midas cichlid fishes in Crater Lake Apoyo an ideal study system to examine the processes leading to ecological specialization and adaptive divergence in sympatry.
Differences in several of these Midas cichlid species' traits are likely adaptations to their divergent ecological niches (e.g., Barlow & Munsey, 1976;Barluenga et al., 2006;Baylis, 1976;Franchini et al., 2014b;Fruciano et al., 2016;Stauffer et al., 2008). Recent studies suggested that divergent body shapes might reflect differences in swimming performances such as cruising swimming and sprinting (acceleration) performance that are particularly important in the limnetic and benthic habitats, respectively, as was shown previously for other fishes (Barluenga & Meyer, 2010;Franchini et al., 2014b;Meyer, 1989, Supporting References, Table 1). Due to its functional importance and significant impact on the daily energy budget, swimming performance is tightly linked to fitness and can play a crucial role in the processes that facilitate sympatric adaptive divergence (reviewed in Higham et al., 2016; further references in Supporting References). In fact, many traits (e.g., fish morphology and physiology) affecting swimming performance underlie developmental constraints and functional trait-offs that allow phenotypes to evolve only within a restricted space along adaptive axes (e.g., either maximizing cruising or sprinting performance, Table 1), promoting rapid adaptive phenotypic divergence (e.g., Langerhans, 2009;Langerhans et al., 2007;Rouleau et al., 2010;Webb, 1984a; further references in Supporting References).
Studies aiming to explore the genotype-phenotype-performance-fitness landscape of incipient ecological speciation can benefit from being conducted in an integrated comparative and functional framework. These can help to identify the ecological, functional and developmental interactions and trade-offs possibly affecting evolutionary trajectories during divergence.

Accepted Article
This article is protected by copyright. All rights reserved In Nicaraguan cichlid fishes, the functional relationship between body shape and swimming performance has not been explored so far, that is, if and how divergent external morphology affects performance and thus probably fitness in different habitats. Here, we aimed to experimentally test whether the benthic A. astorquii and limnetic A. zaliosus species from Crater Lake Apoyo perform differently in terms of cruising and sprinting swimming performance, and if any measured variation in swimming ability is correlated with morphological variation. We also investigated another physiological factor potentially linked to swimming performance and body shape: metabolic rates, i.e., energetic expenditure at rest and during activity, which can be estimated using oxygen consumption (e.g., McDonnell & Chapman, 2016;Peake & Farrell, 2006;Reidy et al., 1995). At the molecular level, we characterized the transcriptomes for gill and muscle tissues of both species at rest and during exercise to explore transcriptional differences between species and activity states and identify the potential underlying key genes and pathways.
According to their body features and habitat use, the benthic high-bodied A. astorquii would be expected to show enhanced sprinting swimming ability, while increased cruising swimming performance is predicted for the limnetic elongated species A. zaliosus. Differences in swimming performance related to divergent habitat usage may further have coevolved with overall metabolic rate: pelagic predators (such as A. zaliosus) often evolved increased metabolic rates and thus energy demand to support cruising swimming. Alternatively, the more elongated morphology of A. zaliosus may reduce drag during swimming and thus also energy demand compared to the deep-bodied A. astorquii. We expect to find these metabolic differences reflected in divergent oxygen demand and gene expression patterns between species. Finally, further species-specific trade-offs may affect these species' metabolic rates, which are currently unexplored in Midas cichlids (e.g., Ellerby & Gerry, 2011;Killen et al., 2016;Webb, 1984b; further references in Supporting References).
This integrated study sheds light on the importance of body shape as an ecologically and evolutionarily relevant trait during niche diversification. We suggest that environmental oxygen availability and heterogeneous habitat use could be additional environmental/physiological variables that play a role in ecological specialization and thereby adaptive divergence in sympatry.

Accepted Article
This article is protected by copyright. All rights reserved We performed four experiments to test for species' differences in swimming performances, metabolic rates and gene expression profiles at rest and during activity. Specifically, the first experiment aimed to investigate swimming abilities and oxygen consumption in cruising swimming (Exp1); the second experiment tested swimming performances during sprinting (Exp2); the third experiment estimated oxygen consumption and gene expression at rest (Exp3); finally, the fourth experiment explored transcriptome divergence during activity (Exp4).

Experimental fishes
One brood from each species of Amphilophus astorquii and A. zaliosus (Fig. 1A,B) were obtained from each one mating pair of laboratory-raised fish collected as fry in Lake Apoyo (Nicaragua) in 2007 with permission of the local authorities (the Ministerio del Ambiente y los Recursos Naturales of Nicaragua, MARENA). Parents were selected randomly to avoid bias and our morphological measurements (see below) confirm that sampled fish reflect the expected natural divergence in morphology between species, as plasticity was found to not affect these fishes' morphology significantly (Franchini et al., 2014b;Kautt et al., 2016b).
The two broods were raised simultaneously in two tanks under the same laboratory conditions. Fish were individually labeled using PIT tags and raised for two years (i.e., until they were all sexually mature). Individuals of the two species were tested in an alternating fashion and all in the same season. To improve concision, more detailed protocols for some analysis steps are available in the Supporting Protocol and Supporting Statistics 1&2. All calculations and statistical analyses were conducted in R v. 3.4.2 (Levin et al., 2016).

Swimming performances and oxygen consumption: protocols in Experiment 1 (Exp1) and 2 (Exp2)
Swimming trials were conducted after a first set of photos for morphometric analyses were taken For each trial, a single fish was tested individually (total n=34, Fig 1D). Body width (BW), wet body mass (BM) and initial measurements of body standard length (iSL) and body height (iBH) were recorded (Fig. 1C). After a preparatory and an acclimation phase, we tested critical swimming speed (U crit ; the maximum speed a fish can swim for at least 20 minutes) or sprint Accepted Article swimming speed (U sprint ; the maximum speed a fish can reach during constant and relatively rapid acceleration of water flow) using two different protocols. These protocols are two of the most widely used methods to measure the cruising and sprinting performance of aquatic organisms, respectively. Briefly, to determine U crit , the water velocity was incrementally increased by 10 cm/s every 20 min until exhaustion (i.e., the fish was pushed against the grid at the end of the swimming chamber) and U crit was calculated according to Brett (1964): U crit = U 1 + [U 2 * (T 1 /T 2 )] in cm/s, where U 1 is the highest velocity maintained for a full 20 min interval (in cm/s), U 2 are the velocity increments between steps (10 cm/s), T 1 is the time spent at the fatigue velocity (in s), and T 2 is the time between increments (1200s). Sprint swimming speed was determined by an immediate velocity increase after acclimatization from 10 to 40 cm/s and subsequent incremental increases by 5 cm/s every 30 s until exhaustion, and U sprint was calculated in the same manner as U crit . Both measurements were corrected for the solid blocking effect using iSL and iBH (more details and references in the Supporting Protocol).
In the first experiment (Exp1), oxygen consumption per kilogram (MO 2 ) was calculated for each fish and velocity phase as the decrease of O2 concentration during the phase and corrected for background respiration (more details in the Supporting Protocol). Despite the fact that oxygen consumption for experiment Exp2 was measured, this data is not presented as the short duration of the protocol may not have allowed the fish to metabolically switch to aerobic muscle movement, which best reflects metabolic rates.
Each fish was tested once in these two experiments (Exp1 and Exp2), allowing a recovery period of at least 24 hours between the two trials. Two fishes (one fish per species) were excluded in Exp1 and Exp2 from downstream analyses of U sprint , U crit and oxygen consumption, as performance results were clearly outliers and swimming behavior was atypical, leading to a sample size of n=16 per species at this point.

Rest/exercise: protocols in experiments Exp3 and Exp4
After a recovery period of several weeks, the remaining fish were split into two groups per species. Only seven individuals per group and species were tested as some PIT tags were lost during the recovery period.
The first group of seven individuals per species was used to estimate oxygen consumption

Accepted Article
This article is protected by copyright. All rights reserved and gene expression at rest (third experiment, Exp3 hereafter, Fig. 1D). Only a very low (5 cm/s) but constant water velocity was used to ensure oxygen replenishment without forcing the fish to swim. Measurements were conducted overnight by alternating phases of oxygen measurement with phases for oxygen replenishment, with all phases taking one hour except a final 1.5h phase oxygen replenishment. Oxygen consumption per kilogram and phase was calculated as in Exp1.
The second group of seven individuals per species was used to obtain exercised fish in oxygenated waters for transcriptome analyses (fourth experiment, Exp4 hereafter, Fig. 1D). For this purpose, water velocity was incrementally increased by 10cm/s every 5min to a speed that all fish were able to swim (according to Exp1; 70 cm/s), which took 35 minutes in total -well beyond the duration for which a switch from anaerobic muscle movement to aerobic muscle movement is expected to occur (Hammer, 1995). The water was also kept at high dissolved oxygen levels (>95%) throughout the experiment (more details on Exp3 and Exp4 in the Supporting Protocol).
After experiments Exp3 or Exp4, fish were immediately euthanized (MS222, 0.4%, pH adjusted). The whole second gill arch and an epaxial white muscle tissue sample were immediately dissected for transcriptome analyses, moved to RNAlater (Merck, R0901) and stored at -20°C.

RNA extraction, sequencing, transcriptome analyses and GO analysis
RNA was extracted from muscle and gill samples, RNA-seq libraries were synthesized and sequenced (151bp, paired end on an Illumina HiSeq X Ten platform). For a few samples, no highquality RNA or cDNA library could be obtained from stored tissue and the total sample size was adjusted to maintain a balanced design with a total of 40 samples (5 sample x 2 species x 2 tissues x 2 activity states; for a detailed protocol of RNA extraction and library synthesis see Supporting Protocol). A reference-based transcriptome was assembled and annotated using the Midas cichlid genome as reference (unpublished data) and transcript per million (TPM) expression tables were generated for each sample (see Supporting Protocol for details). Pairwise analyses of differential gene expression were conducted for each tissue type separately and species-independent transcriptional change between rest and exercise state was calculated using "species" as nested factor. The STRING software (v 11.0) was used to explore putative protein regulatory networks

Accepted Article
This article is protected by copyright. All rights reserved

Morphological data collection, pre-processing and repeatability
Standardized photographs of each tested fish were collected after experiments Exp1 and Exp2 (three photos per individual per experiment). These photos were used to manually determine the coordinates of 17 landmarks (LM), 11 semi (sliding) landmarks and two points that were used to infer the fish's eye centers (Fig. 1 C). Distances were translated to metric units and standard length (SL) was determined as the distance between LM1 and LM6. To assess the repeatability of the morphometric dataset, landmarks were digitalized three times from each photo. Using the gpagen() and procD. Since the obtained morphological variables were initially highly auto-correlated and often confounded with fish size or between species, a number of corrections were applied. All morphological data were corrected for size differences among individuals to remove variation in body shape linked to age (and thus size; Meyer, 1990). As an estimator of size, we initially considered standard length (SL), centroid size and the natural log of individuals' body mass (log(BM)). Overall, conclusions using these different size estimators were similar. However, in our opinion, log(BM) reflects overall size more appropriately and is least affected by body shape differences; thus, for all analyses provided in this manuscript we always refer to the log(BM) results (the full analysis using the more traditional standard length as size estimator is available in the Supporting Statistics 2). To correct for species differences in log (BM)

Accepted Article
This article is protected by copyright. All rights reserved The procD.lm() function (3,000 permutations; R package geomorph) was used to perform Procrustes ANOVAs with permutation procedure. In each of these analyses, we used one experiment-wise coordinates dataset as dependent variables and species, log(BM) and their interaction (Species*log(BM)) as independent variables. Independent variables significantly affected body shape (see Results section for details) in both experiment-wise datasets and thus the residuals of these models were used to perform Principle Component Analyses (PCA).
Subsequently, for each experiment, the corresponding principle components (PCs) individually explaining at least ~5% of total variation were selected, which corresponded to PC1-5 in both geometric morphometric PCAs (Exp1 and Exp2).
To correct body width (BW, a linear morphometric measurement) for the confounding effects of fish size (log(BM)) and species, a linear regression model was fitted for each experiment (Exp1 & Exp2) separately using size (log(BM)) and species (Species) as well as their interaction (log(BM)*Species) as independent variables and BW as dependent variable. As only size affected BW significantly (p=0.011 and p=0.011 for Exp1 and Exp2, respectively, see Supporting Data), the Species term and the interaction terms were removed, and the models recalculated. Downstream analyses used the BW residuals (BW_residuals) derived from these models. Following these steps, our final models were calculated using predictors with very low/no redundancy.

Model selection approach
To explore the relationship between swimming performance (Exp1 & Exp2) or oxygen consumption (only Exp1) and explanatory variables (described in the following paragraphs), a standardized model selection approached using the R stepAIC() function was used to select a suitable linear model based on BIC scores (R package MASS v. 7.3-45;Venables & Ripley, 2002).
The species term, being the focus of this study, was always retained in the final model and an ANOVA type II table was produced using the Anova() function (R package car v. 2.1-2; more details on the model selection approach can be found in the Supporting Protocol).

Effects of the species term and morphological parameters on swimming performance and oxygen consumption
Swimming performance data were parametric in both experiment Exp1 (U crit ) and Exp2 (U sprint ) as

Accepted Article
This article is protected by copyright. All rights reserved suggested by the Shapiro-Wilk test (Shapiro & Wilk, 1965) and the Levene test (Levene, 1960, R package car). The relationship between morphology (selected PCs and BW_residuals), Species and log(BM) (as independent variables), and swimming performance (U crit or U sprint , dependent variable) was investigated separately for the two experiments using univariate linear models for Exp1 and Exp2, respectively. The model selection approach described above was applied to select the most suitable models describing U crit and U sprint (Supporting Data).
After confirming normality of oxygen consumption using the function qqp() with a normal distribution (R package car), we tested the relationship between MO 2-as dependent variable and the same independent variable used for swimming performances tests. As MO 2 was measured during different velocity phases, velocity was included as another independent covariate. To consider the non-independence between the velocity phases, we included fish ID as a random effect. We used linear mixed-effects models fitted by maximum likelihood (R package nlme v. 3.1-128). A similar linear mixed-effects model was used to explore the relationship between MO 2 at rest (Exp3, as dependent variable, normally distributed), and species and "measure phase" as independent variables. The same model selection approach was used to select the most appropriate explanatory variables (Supporting Data). For all final models the Anova() function (package car) was used to compute type II ANOVA tables.
This study was conducted in accordance with local animal welfare regulations and approved by the local Animal Welfare Board (permission number G-16/42 and G-17/54, Freiburg, Germany). All efforts were made to minimize animal suffering through careful collection, handling, and swimming trials based upon the natural rheotaxic behavior and self-motivation of individuals.

Species and size affect morphological variation
We analyzed morphological differentiation between the two focal species using a geometric morphometric approach (Fig. 1 and Fig. S1, S2, S3). Using aligned landmarks, species differed significantly in body shape (Procrustes ANOVA: Species F=27.98, df=1, p<0.001; Fig. 2, Supporting Data). A. astorquii showed the expected larger body-height/body-length ratio (i.e., a deeper body) and a relatively shorter caudal peduncle. It also showed a much more pronounced nuchal hump

Accepted Article
This article is protected by copyright. All rights reserved (a protuberance on the forehead that is prominent in males, clearly visible in A. astorquii in Fig.   1A) when fishes of similar SL are compared between species (Fig.1, Fig.2, Fig. S2). In addition, size (log(BM)) significantly affected body shape in both species (Procrustes ANOVA: log(BM) F=2.36, df=1, p<0.001, Fig. 2, Supporting Statistics 1): with increasing body mass, fish become more deepbodied and the relative length of the caudal peduncle shorter. However, this trend seemed to be less pronounced in A. zaliosus compared to A. astorquii (Procrustes ANOVA: Species*log(BM) F=2.36, df=1, p<0.001; Fig. 2, Supporting Statistics 1, for Exp2 see Fig. S1).

Divergence in cruising but not sprinting swimming performance between species
We investigated the association of species, size, body width and geometric morphometrics data Sprint swimming speed (U sprint , Exp2) was not different between species (F=0.11, df=1, p=0.739; Fig. 2C). In contrast, we identified a significantly (F=15.85, df=1, p<0.001) negative effect of body width on sprint performance, suggesting that wider fish showed lower sprint performance. However, as body width measurements were negatively correlated with PC3 (Supporting Statistics 1), the morphological variation reflected by PC3 may also have contributed to this, suggesting that more anteriorly shifted jaw tips may also be linked to higher sprint performance. A significant negative effect of PC5 (F=5.52, df=1, p=0.03) on sprint swimming speed indicated that fish with longer caudal peduncle showed higher performance.

Higher oxygen consumption in A. zaliosus during activity but not rest
To examine the link between the species term, external morphology and active metabolic rate, we modeled oxygen consumption (MO 2 ) during the critical swimming experiment (Exp1) and at rest (Exp3) using species, the geometric morphometric datasets (PCs), body width and flow velocity as independent variables. Individuals were included as random effects (model selection steps, full statistics and estimates of effects are presented in Supporting Statistics 2). Our model suggests that A. zaliosus consumed slightly less oxygen than A. astorquii at very low velocities (<30 cm/s)

Transcriptomic divergence between species is most pronounced during activity, particularly in gills
The reference-based approach used for our transcriptome analysis allowed us to estimate expression values for the 22,495 genes included in our most recent in-house Midas cichlid genome annotation, although, as expected, not all genes were expressed in each tissue and activity state (Supporting Data). This dataset was used to identify differentially expressed genes (DEGs) between the two tested states (rest vs. exercise in Exp3 and Exp4 respectively) and species (A. astorquii vs. A. zaliosus) for each of the two investigated tissues (gill and muscle). Whole transcriptome tissue-wise PCA analyses suggested that gill samples showed distinct expression patterns according to species and activity state. In contrast, muscle transcriptomes did not ( Between species' gill tissues, 37 genes were differentially expressed at rest (46% of them expressed at higher levels in A. astorquii compared to A. zaliosus) and 68 during exercise (74% showed higher expression levels in A. astorquii than A. zaliosus; Fig. 3A). 18 of these genes were shared between species and 17 of them showed the same direction of expression differences between species at rest and exercise (Fig. 3B). The single shared gene presenting a different direction of regulatory change between states was snf8.

Accepted Article
This article is protected by copyright. All rights reserved We identified DEGs between rest and exercise in Midas cichlid gills by pooling both species gill samples (but accounting for a "Species" effect via a nesting factor). A total of 531 DEGs were identified, 128 of which were upregulated during exercise and 403 up-regulated during rest.
Proteins encoded by the 497 annotated genes among these 531 DEGs showed significant more regulatory interactions than expected by chance (PPI enrichment p-value < 0.001). Additionally, a significant functional enrichment of mitogen-activated protein (MAP) kinase phosphatase genes was detected (FDR p-value = 0.0152). Within species, 146 and 346 DEGs were differentially expressed between activity states for A. astorquii and A. zaliosus, respectively (Fig. 3A). Of these, 12% and 24% showed higher expression at exercise in A. astorquii and A. zaliosus, respectively. 103 DEG were shared between both species DEGs (out of the 146 and 346 DEGs) and 102 of them showed the same direction of regulatory difference (i.e. the direction of their significant regulatory difference between rest and activity does not differ between species; Fig. 3B). Again, one shared gene, snf8, showed the opposite direction of regulatory change from rest to exercise between species.
In muscle tissue, species were differentiated by seven DEGs during rest, of which two showed higher expression levels in A. astorquii compared to A. zaliosus, and eight DEGs during exercise, of which three were more highly expressed in A. astorquii compared to A. zaliosus (Fig.   S9B). Pooling species (but accounting for its effect via a nesting factor) yielded a total of 115 DEGs when comparing rest and exercise, of which 46 were upregulated during exercise and the remaining 69 during rest. 113 of these genes were annotated and STRING analysis again suggested significant more regulatory interactions than expected by chance (PPI enrichment p-value < 0.001). In muscle, a significant functional enrichment of AP1 transcription factor domains was detected (FDR p-value = 0.017). Within species, we found 38 and 23 genes to be differentially expressed between activity states for A. astorquii and A. zaliosus, respectively (Fig. S9A). Of these, 53% and 48% were more highly expressed in resting muscle tissue for A. astorquii and A. zaliosus, respectively. Among the 38 and 23 DEGs between resting and exercise stages in A. astorquii and A. zaliosus, respectively, six were shared between species and all of them showed the same direction of expression differences, i.e. the direction of the significant differences between rest and activity did not differ between species. Three out of the seven and eight genes differentiating the species

Accepted Article
This article is protected by copyright. All rights reserved at rest and exercise, respectively, were shared and the direction of significant difference between species is the same for both activity states (Fig. S9B).
Genes whose expression differed between species exclusively during activity are most likely linked to adaptive differences in swimming physiology, and DEG numbers in gill tissue reflect patterns found for this physiological divergence particularly well (50 genes + snf8, Fig. 3B). The annotation of these genes revealed that many of them were related to gene expression regulation, signal transduction and intra-cellular transport (Supporting Data). These genes might represent basal homeostasis genes. Alternatively, due to the brevity of our exercise protocol (Exp4), these genes might represent the activated early genes of cascades leading to the differential expression of final effector genes, such as genes regulating metabolism or oxygen intake. Specifically, genes that were upregulated in A. astorquii were linked to stimulation of energy expenditure, regulation of feeding and locomotor behavior (such as cart2, Abbott & Volkoff, 2011) as well as molecular responses to increase oxygen availability in tissue such as angiogenesis (tph1a, ece2a, cxc4a, Bussmann et al., 2011;Hyndman & Evans, 2007;Kwan et al., 2016). The one gene that did completely inverse the direction of expression difference between species from rest to exercise, snf8, is likely involved in the regulation of catabolic processes (Malerød et al., 2007;Yeghiayan et al., 1995). Several genes with differential expression patterns between species in gill tissue at either rest (20 genes) or exercise (51 genes) appeared to be associated with differential MO 2 during exercise in the two species (Fig. 3B). To explore co-expression patterns among them, hierarchical cluster analyses were performed on these genes' expression data (Fig. 4). During exercise, eleven genes showed decrease expression in A. astorquii compared to A. zaliosus, while the opposite pattern was found for 39 genes (one gene showed an ambiguous patterns). At rest, opposite patterns were observed with respect to gene numbers: most genes (n=14) had decreased levels of expression in A. astorquii compared to A. zaliosus, while fewer genes showed the reverse pattern (n=5; one gene with an ambiguous expression pattern). Cluster analyses were not performed for muscle due to the low number of identified differentially expressed genes.

Accepted Article
This article is protected by copyright. All rights reserved The Midas cichlid fish benthic-limnetic species pair Amphilophus astorquii and A. zaliosus inhabiting the Nicaraguan Crater Lake Apoyo represents a striking example of adaptive divergence with gene flow along the benthic-limnetic axis (Barluenga et al., 2006;Elmer et al., 2014;Meyer, 1990). We tested these two species for differences in swimming performances that are thought to be adaptive in their respective divergent ecological niches (benthic vs. limnetic; Supporting References). Our results indicate that these species differ in cruising swimming performance, oxygen demand and gene expression during activity. We provide empirical evidence that external body morphology together with physiological mechanisms indeed constitute ecologically relevant traits that likely represent adaptations to the limnetic and benthic habitats. These traits facilitate niche segregation and the inferred performance trade-off likely contributes to adaptive divergence further supporting the hypothesis of ecological speciation.

Divergent morphology and swimming performance in Midas cichlids
Body shape variation between A. astorquii and A. zaliosus has been shown to mirror ecological specialization in habitat use (Barluenga et al., 2006;Franchini et al., 2014b). In line with these previous studies, we confirmed that the two Midas species differ in their external morphology, most pronouncedly in body depth (A. astorquii being more high-bodied but shorter), body mass to body length ratio (A. astorquii being considerably heavier per cm standard length) and the nuchal hump size (see also Methods and Results sections).
Consistent with ecological differences and predictions on functional morphology (Table 1, Supporting References), we found a higher cruising performance in A. zaliosus than in A. astorquii.
This difference may be due to divergent body morphology (e.g. resulting in different drag, Fig. 1), divergent physiologies (e.g. divergent metabolic rates, Fig. 3) or a combination of both. We observed species' differences in nuchal humps, a male sexual trait that is strongly preferred by females in Midas cichlid fishes (Barlow & Munsey, 1976;Barlow & Siri, 1997) but could be costly in terms of swimming performance (the handicap principle; Zahavi, 1975). However, divergent critical swimming speed is established even in individuals without pronounced nuchal humps (i.e., our body mass-matched sub-selection). No other morphological variables were significantly linked with critical swimming performance indicating that, if external morphology affects critical swimming ability, this variation is already accounted for in the species term of our models, which Accepted Article includes other species-related factors such as divergent physiologies. Despite this confounding effect of species on body morphology, it is reasonable to assume that the observed difference in performances is at least partially due to divergent body shapes, as differences in external body shape have been linked to divergent swimming performances before (Barluenga et al., 2006;Franchini et al., 2014b, Supporting References, Table 1). This supports the notion that A. zaliosus's distinct limnetic morphology might have evolved specifically in response to selection for increased cruising performance in their open water habitat, an ability that would be expected to enhance their catching ability of evasive and patchily-distributed food sources (Barluenga et al., 2006;Elmer et al., 2014;Franchini et al., 2014a).
The benthos-associated fish A. astorquii was expected to outperform A. zaliosus in sprinting swimming performance. However, we did not observe significant differences in sprint swimming speed (Exp2) between these two species. Although these cichlid species exhibit conspicuous differences in their external morphologies, swimming abilities are the result of the combination of numerous factors including body shape, physiology, overall condition and swimming behavior. This complex physiological integration could have masked more subtle differences, decoupled physiological systems or the use of different underlying physiological or locomotion strategies (e.g., Crespel et al., 2017;Hendry et al., 2011;Law & Blake, 1996).
From a morphological perspective, sprint swimming performance (Exp2) was negatively correlated with body width (BW, which did not differ between species, see Methods and Results sections), potentially resulting from an increase of drag associated with this feature (Van Wassenbergh et al., 2015). Mechanistic trade-offs may not allow these species to evolve a further reduced body width as it may compromise the functionality of other traits such as pharyngeal jaw crushing or gill oxygen uptake (Chapman et al., 2000). Sprint swimming performance was also positively associated with caudal peduncle length (PC5) in association with almost no other morphological differences (except for slightly more dorsally operculum, LM 13 and 27; PC5). This finding seems to suggest that increased caudal peduncle length may lead to enhanced propulsion.
Interestingly, A. zaliosus typically have relatively longer caudal peduncles, but we did not observe differences between species in sprint swimming test (Exp2).

Accepted Article
This article is protected by copyright. All rights reserved Oxygen consumption (MO 2 ) was used as a measure for energy demand and thus metabolic rate during cruising swimming and rest. We expected divergent oxygen demand between the streamlined limnetic A. zaliosus and the deep-bodied benthic A. astorquii; however, the direction of divergence was hard to predict as opposing effects have been observed in fishes: while streamlined species produce less drag than deep-bodied ones and thus should need less energy (and thus oxygen) to swim, pelagic swimmers tend to generally have higher metabolic rates to support cruising swimming (Ellerby & Gerry, 2011;Killen et al., 2016), which predicts higher oxygen demand. In line with the first hypothesis, more pointy-faced fish with shallower bodies required less oxygen (PC4 of Exp1). However, we also documented increased MO 2 in A. zaliosus compared to A. astorquii in cruising swimming (Exp1) and, independently of species, fish of larger size (i.e. mass) have a lower MO 2 per kg body mass during exercise than smaller individuals, likely reflecting size-dependent relative drag reduction (Hoerner, 1965;Webb, 1975). Alternatively, the relationship of metabolic rate to body weight may in this species deviate from the classically assumed MO 2~B M -1 , which we also applied in this study (see Supporting Protocol). However, the used metabolic scaling is unlikely to affect our species comparisons as we used fishes of similar weight. Future studies may thus strive to determine the precise relationship between body mass and metabolic rate, which would allow studies comparing metabolic rates using fish of different sizes (and thus body mass).
Divergent MO 2 during activity may represent a physiological adaptation to different oxygen concentrations where A. astorquii and A. zaliosus live and swim most actively (feeding). Trophic analyses revealed that prey and foraging grounds differ markedly between these two species (Fig. 5, Barluenga et al., 2006;Elmer et al., 2014;Franchini et al., 2014a). On the one hand, A. zaliosus exhibits trophic signatures typical for the limnetic habitat closer to the water surface that is enriched in oxygen and food of higher quality compared to the benthic environment. On the other hand, A. astorquii mainly feeds in the benthic zone at lower trophic levels (Fig. 5); only its stomach contains Chara algae. These algae are bottom-dwellers that prefer clear, but deeper and less oxygenated waters (Barluenga et al., 2006;Elmer et al., 2014;Guha, 1995). Accordingly, the selective pressures to cope with low oxygen levels while actively foraging might have been more pronounced in A. astorquii than A. zaliosus. Trade-offs between swimming performance and tolerance to low resource availability, particularly oxygen, are known from other teleosts. They

Accepted Article
This article is protected by copyright. All rights reserved contribute to divergent selection pressures that strongly favor ecological segregation: high metabolic rates and increased cruising performances in pelagic species and decreased metabolic costs and performances in benthic species ("cost reduction strategy", Killen et al., 2016). Our findings suggest that this pattern of ecological differentiation might be observed also in this Midas cichlid species pair (Fig. 5).
We observed no divergence in metabolic rate between species at rest, but A. zaliosus consumed significantly more oxygen with increasing velocity. Habitats used for resting and feeding are not only different for the species (Barluenga et al., 2006;Torres-Dowdall, personal communication), but also likely more heterogeneous for A. zaliosus compared to A. astorquii, possibly favoring a more variable metabolism to cope with these environmental differences in the former (Fig. 5). In fact, as the former forage rather in open waters, they might need to descend to considerable depths to rest where dissolved oxygen levels are lower, while rest and feeding occurs at similar depths in the benthic species. Recent environmental measurements in Lake Apoyo showed that with increasing depth the dissolved oxygen levels are relatively constant at ~90% until at 10m depth a sudden drop to ~60% occurs (Torres-Dowdall & Meyer, In press).
Heterogeneity in their habitat may have thereby contributed to the evolution of metabolic plasticity in A. zaliosus.

Differentiation in gill and muscle transcriptomes
To our knowledge, this study is among the firsts to explore divergence in gene expression between recently diverged species and activity states (rest and exercise, Exp3 and Exp4 respectively) in target tissues (but for trout see e.g. Magnoni et al., 2013;Palstra et al., 2013). A PCA of gill and muscle transcriptomes suggests that gene expression was markedly differentiated between species and activity states in gills, and in contrast to muscle that showed no clear differences (Fig. S9). A. astorquii and A. zaliosus are closely related species that diverged less than 24,000 years ago (Barluenga et al., 2006;Kautt et al., 2016a); strong divergent selection is thus expected to have caused these differences in genes expression. Alternatively, the observed pattern might reflect drift affecting tissue-specific expression or compensatory gene regulation.
When all white muscle samples per activity state were pooled (and species used as nested effect), 115 genes were differentially expressed and these were enriched in AP1 transcription factor Accepted Article binding domains. Limited differentiation in muscle tissue within and between species may be due to the usage of white muscle (instead of red muscle) tissue. While white muscle is thought to be responsible for short-term burst swimming and less reliant on high oxygen supply, red muscle allows for endurance swimming at the cost of higher oxygen demand (Videler, 2012). It is thus possible that red muscle transcriptomic patterns would have shown a more pronounced divergence between species than white muscle did. However, a previous study exploring divergence between rest and exercised red and white muscle transcriptomes of trout found overall very similar differentiation magnitudes between muscle types and activity states across genes (Palstra et al., 2013). Thus, both red and white teleost muscle may not respond strongly to the activity state by altering the expression of specific gene groups in Midas cichlids.
Overall, our PCA analysis for gill tissue suggests that the transcriptomic change between the rest and exercise state was the main source of variation in the transcriptomic dataset ( Fig. S8: PC1, differentiating activity states, reflects ~29% of total variation), while species differences came second ( Fig. S8: PC2, differentiating species, reflects ~14% of total variation). The relatively small number of upregulated genes during exercise (n=128) compared to downregulated genes (n=403) suggests that during exercise transcriptional rates were decreased more broadly while genes whose expression increased contributed to more specific functions. The regulatory interaction analysis (Fig. S10) suggests that more interaction between proteins encoded by DEGs were found than expected by chance, which would be expected when few, specific regulatory pathways underlie the transcriptomic difference between states. A prime candidate for such a pathway is the MAP kinase pathway, as it was identified to be significantly enriched among DEGs. MAP kinase signaling is critically involved in signal transmission from the environment into the cell and thus in regulating many cellular processes (Seger & Krebs, 1995) and it has been shown to respond to oxidative stress as well as physical exercise (Aronson et al., 1997;McCubrey et al., 2006).
Specifically, all dual specificity phosphatase genes among DEGs (dusp1, dusp2, dusp5, dusp8a, dusp16) were down-regulated in exercise samples. These genes inactivate MAP kinases by dephosphorylation of either MAP kinase tyrosine or threonine residues (Camps et al., 2000), which suggests that MAP kinase signaling is up-regulated during exercise. Our STRING analysis also revealed that the largest identified regulatory cluster contained numerous genes contributing to cell proliferation or control the cell cycle (Supporting Data): while genes controlling cell arrest

Accepted Article
This article is protected by copyright. All rights reserved were down-regulated during exercise (e.g. gadd45b and gadd45g) those assisting DNA replication were up-regulated (e.g. mcm5, mcm6 and mcm7; Kearsey & Labib, 1998;Marhin et al., 1997). We conclude that gill tissue responded to physical exercise in our experiment likely by increasing gill cell proliferation and differentiation. For muscle tissue similar but less pronounced patterns were identified. A possible major regulatory pathway mitigating the white muscle response to exercise might be the AP1 transcription factor (Fig. S11), which has been described to be involved in mitigating oxidative stress in response to exercise in skeletal muscle tissue (Hollander et al., 2001;Zhou et al., 2001).
At the species level, patterns of species' differentiation in the gill transcriptome suggest that the transcriptomic change from rest to exercise is more pronounced in A. zaliosus than in A.
astorquii and that differentiation between species is more pronounced during exercise than rest, similarly to the pattern observed for MO 2 . Complementarily, genes that are differentially expressed between species during activity involve the regulation of energy expenditure and response to hypoxia. Considering that water was constantly oxygenated during the exercise protocol (Exp4), the observed expression patterns suggest that a genetic adaptation to divergent oxygen levels evolved in response to different oxygen concentrations that are known to exist at different depths in Lake Apoyo. Thus, our data suggests that physiological and molecular adaptations in response to divergent oxygen levels associated with habitat specialization have likely coevolved along with feeding ecology of A. astorquii and A. zaliosus.

Mechanisms underlying adaptive divergence
Mechanisms such as ecologically-based disruptive selection and assortative mating counteract gene flow and may lead to sympatric speciation (Bolnick, 2011;Gavrilets et al., 2007). It is currently not feasible to directly measure the strength of selective forces affecting the benthiclimnetic divergence in Lake Apoyo (as described in Kautt et al., 2016b) or the effect of performance on (ideally, lifetime, see Arnold, 1983) fitness in the field as these species are endemic to an habitat that is not easily accessible for detailed, long-term observations or experiments in large sample sizes. Nonetheless, in line with earlier findings (Barluenga & Meyer, 2010;Barluenga et al., 2006;Elmer et al., 2014;Elmer et al., 2010;Franchini et al., 2014a;Franchini et al., 2014b;Fruciano et al., 2016;Kautt et al., 2016a;Meyer, 1990;Stauffer et al.,

Accepted Article
This article is protected by copyright. All rights reserved 2008), this study provides indirect empirical support for strong ecologically-based disruptive selection, likely associated with competition and performance trade-offs leading to habitat specialization in Lake Apoyo (Fig. 5).
In conclusion, ecologically-based divergent selection on swimming performance and oxygen metabolism and trade-offs in key ecological traits, together with other mechanisms such as genomic linkage, habitat preferences and assortative mating (Barluenga et al., 2006;Baylis, 1976;Elmer et al., 2014;Elmer et al., 2010;Franchini et al., 2014b;Fruciano et al., 2016;Kautt et al., 2016b;Stauffer et al., 2008), indeed likely promoted differentiation in this sympatric Midas cichlid species pair in Lake Apoyo. Future studies replicating our integrated experimental design in other families of the same species pair or sympatric species pairs, such as A. amarillo and A. sagittae in the Nicaraguan crater Lake Xiloá, will further clarify the processes enabling divergence despite gene flow in these cichlid fishes as well as other iconic occurrences of sympatric ecological speciation.
width, wet mass), geometric morphometric coordinates (raw and processed), swimming performance data for each individual and oxygen consumptions during the critical swimming

Author Contributions
FR, RFS, PF and AM designed the study. FR and AM obtained the necessary grants, equipment and permissions. FR assembled the swimming tunnel system and optimized swimming performance protocols. Swimming performance data was collected by FR and analyzed by FR and RFS.
Morphological data was collected by FR and RFS and analyzed by RFS and FR. Transcriptome data was collected and analyzed by RFS and PF. AFK provided the annotation for the transcriptomes. FR and RFS drafted the manuscript. All authors revised and agreed to the manuscript. Table 1 Investigated swimming performances, expected associated morphological characteristics enhancing these abilities, predicted habitat-linked effects on fish fitness and protocols used to measure these swimming performances in this study (Higham et al., 2016;Langerhans, 2009;Langerhans et al., 2007;Rouleau et al., 2010;Webb, 1984a