Repeated divergence in opsin gene expression mirrors photic habitat changes in rapidly evolving crater lake cichlid fishes

Selection pressures differ along environmental gradients and organisms’ phenotypes. Traits tightly linked to fitness (e.g., the visual system) are expected to closely track environmental variation along gradients. Within such gradients, adaptation to local conditions might be due to heritable and non-heritable, environmentally induced variation. Disentangling these sources of phenotypic variation requires studying, in nature and the laboratory, closely related populations experiencing different environments. The Nicaraguan great and crater lakes show an environmental gradient in photic conditions extending from clear crater lakes to very turbid great lakes. From two old, turbid great lakes, Midas cichlid fish (Amphilophus cf. citrinellus) independently colonized seven isolated crater lakes of varying light conditions, resulting in a small adaptive radiation. We estimated the variation in visual sensitivities along this photic gradient by measuring differential cone opsin gene expression among populations from different lakes. The visual sensitivities observed in all seven derived crater lake populations have not changed randomly but shifted predictably in direction and magnitude, repeatedly mirroring changes in photic conditions. Intrapopulation phenotypic variation decreases as environments become spectrally narrower suggesting different selective landscapes within the gradient. Comparing wild-caught and lab-reared fish revealed that 48% of this phenotypic variation is genetically determined and evolved rapidly. Our results demonstrate deterministic, rapid phenotypic evolution that fine-tunes visual sensitivity to fine-scale environmental variation.

Methods
Study Design
In Nicaragua, a natural experiment occurred, where from a common source population in the great lakes Managua and Nicaragua, seven isolated crater lakes were independently colonized by Midas cichlid fish (Amphilophus cf. citrinellus) between 4,700 to 800 years ago (Barluenga et al. 2006; Kautt et al. 2016; 2020). The young radiation of Midas cichlids currently encompasses 13 nominal species characterized by genomic and morphological differentiation among lakes with sympatric and allopatric species showing divergence in traits related to lip size, body shape, pharyngeal morphology, or body coloration (Fig 1a, Torres-Dowdall and Meyer 2021). Given the isolated nature of the crater lakes, their colonization from common source populations and their geomorphological similarities, crater lakes have been considered as natural replicates regarding many of their environmental factors (Kautt et al. 2018). However, photic conditions differ widely among Nicaraguan lakes. The great lakes are big and shallow and the winds create waves that constantly stir up the sediments making them very turbid (Elmer et al. 2010). Similar light conditions can be found in the San Juan River, a major river connected to Lake Nicaragua that is also inhabited by Midas cichlids. In contrast, the crater lakes are very deep and thus sediments are deposited at depths far from the influence of waves which might contribute to the observation that most crater lakes are clearer than the great lakes, but there is still substantial variation among crater lakes (Torres-Dowdall and Meyer 2021). Studies on the visual ecology of the two oldest crater lakes Apoyo and Xiloá found convergent changes in the visual system of these populations regardless of their color morphs suggesting a strong effect of the ambient light environment (Härer et al. 2018; Torres-Dowdall et al. 2017). However, given the diversity of photic environments found among crater lakes, their population isolation (i.e., absence of gene flow), and their recent and independent colonization from a common source, the system provides an excellent opportunity to investigate the predictability of phenotypic evolution in visual systems along a wide range of photic conditions by asking, how closely do visual systems track local conditions along a natural gradient?

Characterization of the Photic Environments
To characterize the photic conditions found among the Nicaraguan lakes, underwater irradiance was measured from both great lakes, Managua and Nicaragua, seven neighboring crater lakes, and one riverine population (Fig. 1b) Absolute irradiance was measured using a spectrometer (FLAME-S-XR1-ES, Ocean Insight, USA) connected to a 25-m UV-VIS optical fiber (OCF-104472, Ocean Insight, USA) with a cosine corrector (CC-3-UV-S, Ocean Insight, USA). Multiple consecutive measurements during daytime between 10 am-2 pm were performed at 0.15 m, 1 m, 3 m, 5 m, 10 m, 15 m, 20 m, and 25 m depth. Sites with less than 25 m depth were measured until their deepest point. The measurements were performed by orienting the sensor upwards (downwelling irradiance, Ed), sideways (sidewelling irradiance, Es), and downwards (upwelling irradiance, Eu). Absolute irradiance measurements were corrected for integration time and converted to photons/cm2/s/nm (E) based on Johnsen (2012). To minimize the effect of outliers due to handling of the spectrometer, median absolute irradiance of 3-10 measurements for each depth was used and smoothed using a rolling mean over 5 nm following the manufacturer’s instructions (oceaninsight.com). Only wavelengths within the visible spectrum (350-700 nm) were used based on the peak sensitivity of visual pigments in fish (Carleton et al. 2020; Rennison et al. 2016). To allow the comparison of spectral shape across sites, absolute spectra were divided by their respective maximal value resulting in normalized irradiance.

Opsin Gene Expression
To determine the degree of variation in cone opsin expression in populations of Midas cichlids, six to eight wild-caught adult fish per site were collected in January and February 2018 from ten locations across Nicaragua for a total of 78 individuals (Table S1). In the rest of the text, the term “population” is used to refer to this level of sampling (i.e., locations), although some of the sampled populations correspond to formally described species within the Midas cichlid radiation. Additionally, 62 lab-reared adults from nine populations raised for at least two generations in the animal research facility at the University of Konstanz were included in our study. These laboratory experiments were done to measure phenotypic variability in the absence of developmental noise due to the different light conditions fish experience in the wild. Only sexually mature fish (at least 2y old) were used given that opsin gene expression varies during ontogeny eventually reaching a developmental plateau at adulthood (Härer et al. 2017; 2019). Fish were sampled during the same daytime (11-15 pm) to control for diurnal variation in gene expression (Yourick et al. 2019), and euthanized by applying an overdose of MS-222 and subsequent cervical dislocation. The retinas were removed and stored in RNAlater (Sigma-Aldrich, USA) at -20°C until extraction. RNA was extracted using a standard Trizol-chloroform protocol based on Rio et al. (2010). For each sample, 200 ng of total RNA were used to synthesize first-strand cDNA using the manufacturer’s protocol (GoScriptTM Reverse Transcription System, Promega, USA). Gene expression of six cone opsin genes (sws1, sws2b, sws2a, rh2ab, rh2b and lws) and two reference genes (gapdh and imp2) was measured using quantitative real-time PCR (qPCR) for 40 cycles (CFX96TM, Bio-Rad Laboratories, USA) following Härer et al. (2017). Expression of the paralog rh2aa was not measured since it is not expressed in Midas cichlids (Torres-Dowdall et al. 2017). Mean threshold cycle (Ct) values from three technical replicates were used for analysis. Primer sequences, amplification efficiencies, and mean expression of reference genes are reported in the supplementary material (Table S2-S3, Fig S3). Proportional opsin expression for each individual was calculated as the amount of each cone opsin (Ti) relative to the total cone opsin expression (Tall) as Fuller et al. (2004).

Testing the Association between Photic Changes and Sensitivity Shifts
To determine if changes in the visual system between ancestral and derived populations are correlated with differences in the photic conditions between the source and the derived habitats, we took advantage of our knowledge about the demographic history of Midas cichlids (Kautt et al. 2020). Extensive population genomic analyses showed that Crater Lake Apoyo was colonized from Great Lake Nicaragua and all other crater lakes from Great Lake Managua. While the ancestry of Crater Lake Masaya is admixed, Great Lake Managua is considered its main population source for analysis in this study due to the small contribution of Great Lake Nicaragua (~22%, Kautt et al. 2018; 2020). We asked if photic habitat differences drive the phenotypic divergence seen between source and derived populations of Midas cichlids. For this, we calculated the degree of correlation between spectral attenuation coefficients (Kd) and estimated spectral sensitivity curves (∆SSCj). Spectral attenuation coefficients represent the extinction of ambient downwelling light with depth and are less prone to background noise (e.g., atmospheric events, waves) than irradiance, making it a more robust estimate of the spectral characteristics of water bodies (Mobley 1994; Rennison et al. 2016). The localized spectral attenuation coefficient, Kd, was calculated based on Sabbah et al. (2011). Next, we determined the change in predicted visual sensitivity experienced by fish from derived populations compared to the source populations. For this, the individual continuous estimates of visual sensitivity, the spectral sensitivity curves (SSCj) were estimated for each specimen following Rennison et al. (2016). In short, absorbance templates from Govardovskii et al. (2000) and absorption peaks from Torres-Dowdall et al. (2017) for each opsin were used. We simulated scenarios assuming either only A1- or only A2- chromophore usage. The sensitivity curve was then weighted by the proportional expression of each opsin and the sensitivity curves of the six expressed opsins were added (Fig.S5). Subsequently, shifts in spectral sensitivity curves (∆SSCj) were estimated as the difference between the median spectral sensitivity of the source population and the individual spectral sensitivities of each fish in the derived populations as in Rennison et al. (2016). Finally, we ran Pearson’s correlation tests between shifts in spectral sensitivity (∆SSCj) and changes in attenuation coefficients (∆Kd) for each individual. 

Predictability of Visual Sensitivity based on Photic Conditions
To test if the spectral sensitivity of fish can be predicted based on their local photic conditions, we used point estimates for visual sensitivity and regressed on a composite axis of the photic conditions at their lake of origin at one meter depth. One meter depth was chosen as a compromise to include all sites (shallow habitats ~2m depth) and in agreement with the habitat ecology of Midas cichlids (Dittmann et al. 2012; Oldfield et al. 2006). The composite axis was generated using correlation-based principal components analysis (PCA) of seven z-standardized variables. We used down- and sidewelling λP50 as the spectrum-halving wavelength which summarizes the photon distribution into a single value indicating short- or long-wavelength predominant spectra (McFarland and Munz 1975a). We also included λP25 and. λP75, the wavelengths within which 50% of the photons are found (i.e., spectral broadness). Finally, the percentage of downwelling photons available at one meter compared to 15 cm below the water surface, an estimate of luminosity %Ed, was also included in the PCA (Table S4). The response variable predicted sensitivity index (PSI) was defined as the sum of peaks in absorption of each opsin weighted by its proportional expression in the retina (Hofmann et al. 2009) and calculated using the following equation. To determine the genetic component of the phenotypic variation found in the visual system of Midas cichlids, we also estimated the predicted sensitivity index of individuals reared in the laboratory. Given that phenotypic divergence in the wild could also be mediated by environmentally induced changes, phenotypic variation measured under common-garden conditions informs about its genetic component. We used a linear mixed-effect model to determine the percentage of phenotypic variation that is explained by native photic conditions independently of rearing conditions. The model considered the predicted sensitivity index (PSI) as the response variable, photic environment (PC1), rearing environment (i.e., wild or lab), and their interaction as predictor variables (fixed effects); lake of origin was used as a random intercept. The relative importance of each regressor to the amount of explained variance was estimated based on Stoffel et al. (2021). Confidence intervals for mean regression lines accounting for standard error of the regression line and intercept on each population were calculated following Breheny and Burchett (2017). To account for potential bias in predicting spectral sensitivity from gene expression data, a sensitivity analysis to assess robustness of predictors and outcome to unobserved confounding factors was performed as in Cinelli and Hazlett (2020). Additionally, we regressed the coefficients of variation in our estimates of visual sensitivity (PSI) within each population against the photic axis (PC1) to test if phenotypic variation is reduced in spectrally narrow environments. Diagnostic plots are provided in the supplemental materials (Fig. S10-11).