Multispecies colour polymorphisms associated with contrasting microhabitats in two Mediterranean wrasse radiations.

Intraspecific colour polymorphisms (CPs) present unique opportunities to study fundamental evolutionary questions, such as the link between ecology and phenotype, mechanisms maintaining genetic diversity and their putative role in speciation. Wrasses are highly diverse in ecology and morphology and harbour a variety of colour-polymorphic species. In the Mediterranean Sea, wrasses of the tribe Labrini evolved two species radiations each harbouring several species with a brown and a green morph. The colour morphs occur in complete sympatry in mosaic habitats with rocky outcrops and Neptune grass patches. Morph-specific differences had not been characterized yet and the evolutionary forces maintaining them remained unknown. With genome-wide data for almost all Labrini species, we show that species with CPs are distributed across the phylogeny, but show evidence of hybridization. This suggests that the colour morphs are either ancient and have been lost repeatedly, that they have evolved repeatedly or have been shared via hybridization. Focusing on two polymorphic species, we find that each colour morph is more common in the microhabitat providing the best colour match and that the morphs exhibit additional behavioural and morphological differences further improving crypsis in their respective microhabitats. We find little evidence for genetic differentiation between the morphs in either species. Therefore, we propose that these colour morphs represent a multi-niche polymorphism as an adaptation to the highly heterogeneous habitat. Our study highlights how colour polymorphism (CP) can be advantageous in mosaic habitats and that Mediterranean wrasses are an ideal system to study trans-species polymorphisms, i.e. polymorphisms maintained across several species, in adaptive radiations.

INTRODUCTION

Colour polymorphisms (CPs), the occurrence of two or more distinct, genetically determined colour morphs within a single population (Huxley, 1955), are widespread in nature (Gray & McKinnon, 2006; Limeri & Morehouse, 2014; Stevens, 2016) and contribute a substantial amount of phenotypic variation (Limeri & Morehouse, 2014; McLean & Stuart‐Fox, 2014). CPs present unique opportunities to study the evolution and maintenance of intraspecific variation (Gray & McKinnon, 2006; Jamie & Meier, 2020; Svensson, 2017). A CP may serve predator avoidance (Karpestam et al., 2016), thermoregulation (Rao & Mendoza‐Cuenca, 2016), or communication (Endler, 1978; Stuart‐Fox & Moussalli, 2009), and may allow polymorphic species to exploit a broader niche, whereby distinct colour morphs occupy alternative fitness peaks (Levene, 1953; Maynard Smith, 1962; Svensson, 2017). A genetically determined stable CP requires some form of balancing selection to maintain equivalent fitness among colour morphs (Arias et al., 2016; Gray & McKinnon, 2006; Maynard Smith, 1962; Pfennig et al., 2010; Svensson, 2017), but identifying the ecological or sexual factors exerting such balancing selection is often challenging (Bolton et al., 2016; Forsman, 2016).

Different forms of balancing selection can maintain a CP such as negative frequency‐dependent selection, heterozygote advantage or fitness variation in time and space (Gray & McKinnon, 2006; Jamie & Meier, 2020; Llaurens et al., 2017). Negative frequency‐dependent natural selection may be driven by predator avoidance or competition for mates or territories (Gray & McKinnon, 2006). CP can reduce visual predation if predators form search images for the most common morph, leading to a fitness advantage of rare morphs (Gray & McKinnon, 2006; Karpestam et al., 2016; Maynard Smith, 1962; Sinervo & Svensson, 2002). Sexual selection may also be negative frequency‐dependent with females preferring rare male morphs (O'Donald & Majerus, 1988), males being more aggressive to common male morphs or their own morph (Seehausen & Schluter, 2004), or polymorphisms linked to alternative mating strategies (Sinervo & Lively, 1996). Such CPs maintained by sexual selection or sexual conflict are typically limited to one sex (Ajuria Ibarra & Reader, 2014; Kingston et al., 2003; Outomuro et al., 2014). Alternatively, CPs can be maintained due to heterozygote advantage (Gemmell & Slate, 2006; Kellenberger et al., 2019) or if contrasting microhabitats favour different morphs (multi‐niche polymorphisms, Gray & McKinnon, 2006; Maynard Smith, 1962; Rueffler et al., 2006; Skulason & Smith, 1995). In such multi‐niche polymorphisms, balancing selection arises from competition within each niche, leading to increased fitness of a morph that is rarer relative to the carrying capacity of the niche in which it is favoured.

Distinct colour morphs are not necessarily genetically determined, but may be induced by environmental cues or depend on the developmental stage of an organism (so‐called polyphenisms, Mayr, 1963; Pfenning et al., 2010). In such cases, there is no balancing selection. However, as plasticity is expected to be costly, the alleles underlying a polyphenism should only persist if the ability to express alternative phenotypes is advantageous (Scheiner, 1993).

In many cases, colour morphs also differ in other traits such as life history, behavioural or morphological traits (Gray & McKinnon, 2006; Svensson, 2017). A nearly perfect match between colour morphs and associated traits may exist due to plastic responses of some or all traits to the same environmental cues, pleiotropic effects of a single locus, or a supergene with multiple tightly linked genes (Gray & McKinnon, 2006; Jamie & Meier, 2020; Svensson, 2017). However, if the different morph‐associated traits are determined by physically unlinked genes, matings between morphs produce offspring with suboptimal trait combinations, which may impose selection for morph‐assortative mating potentially leading to speciation (Gray & McKinnon, 2006; Losos & Mahler, 2010; McKinnon & Pierotti, 2010).

Polymorphisms are often associated with species radiations whereby multiple closely related species harbour similar intraspecific morphs (Enbody et al., 2021; Gray & McKinnon, 2006; Hugall & Stuart‐Fox, 2012; Jamie & Meier, 2020). Independent maintenance of the same CP in multiple species suggests a fitness advantage for the rarer morphs or heterozygotes. An association between intraspecific polymorphisms and high speciation rates can arise if morphs act as precursors to speciation (Gray & McKinnon, 2006; Hugall & Stuart‐Fox, 2012; Jonsson, 2001; Seehausen & van Alphen, 1999). Alternatively, the association can be due to the presence of closely related taxa in a species radiation, allowing the (re‐)introduction of advantageous alleles underlying the morphs via gene flow (Jamie & Meier, 2020). In order to understand the relationship between intraspecific polymorphisms and adaptive radiation in a specific study system, it is important to test if the polymorphism represents an ancestral condition, if it evolved independently in several species or if it has been transferred between species by introgression. This requires a robust phylogeny, information of the presence or absence of a CP in each species and tests for introgression.

Wrasses (family Labridae) provide an ideal opportunity to investigate the complex interaction of ultimate drivers and proximate mechanisms maintaining a CP. In the Mediterranean Sea, wrasses of the tribe Labrini have undergone two species radiations in the genera Symphodus and Labrus (Hanel et al., 2002), with several species in each radiation displaying two distinct morphs of either brown or green colouration (Louisy et al., 2015). All polymorphic wrasses are found in the Mediterranean Sea, whereas Atlantic wrasses do not harbour colour morphs. The colour morphs may thus have evolved as an adaptation to the habitat specific to the Mediterranean Sea. Due to the absence of strong tidal action, the shallow parts of the Mediterranean Sea are covered by mosaics of Neptune grass (Posidonia oceanica) patches and rocky substrate. Little is known about the selective forces or the genetic mechanisms underlying the colour morphs in Mediterranean wrasses (Zupo & Stuebing, 2010).

Recent research focusing on the ballan wrasse (Labrus bergylta) revealed genome‐wide differentiation between sympatric plain and spotted morphs in some parts of the distribution range, consistent with early stage speciation (Casas et al., 2021; Quintela et al., 2016). Although no study has investigated genomic differences between green and brown morphs in wrasses so far, aquarium experiments by Arigoni et al. (2002) showed that the CP in Symphodus ocellatus and Symphodus roissali is plastic to some extent with colour changes within several weeks to months, but not plastic in Symphodus rostratus. Intermediately coloured individuals exist within several of these polymorphic species, but it is unclear whether they represent a transitional colour stage (as a result of plasticity) or if they are in fact genetically intermediate.

Here, we first assess the phylogenomic relationships of Labrini species. Next, focusing on the two most common Symphodus species, S. roissali and S. rostratus, we ask whether green and brown colour morphs differ in traits other than colouration and what evolutionary mechanisms are likely to maintain the polymorphism. We find that the colour‐polymorphic Symphodus species are not each other's closest relatives indicating repeated origins or losses of the CP, whereas in Labrus, the three polymorphic species form a monophyletic clade with possibly ancestral CP. In both focal species, we find that the green morph is tightly associated with P. oceanica, and the brown one with rocks covered in brown algae (Cystoseira spp.) forests. Furthermore, the morphs also differ in morphology, aggression and anti‐predator behaviours, all of which are consistent with adaptation to contrasting selection pressures in the two microhabitats. We find no association of colouration with sex or age in either species, and morphs show little, if any, genome‐wide differentiation. Our findings are thus consistent with maintenance of the CP as a multi‐niche polymorphism that is an adaptation to the highly heterogeneous mosaic habitat of rocky shores in the Mediterranean Sea. The many colour‐polymorphic species in both Labrini radiations are thus a great study system to investigate the role of polymorphisms in adaptive radiations.

METHODS

Phylogenetic tree based on genome‐wide data

RAD sequencing

To infer the phylogenetic relationships of Labrini species, we obtained 101 samples covering 21 of the 24 known species (see Tables S1 and S2). For our focal species, we collected a total of 41 S. rostratus and 12 S. roissali samples for population genomic analyses (Table S3). DNA was extracted from fin clips or muscle tissues stored in 100% ethanol or RNA later using a phenol chloroform protocol (Sambrook & Russell, 2001). We prepared six restriction‐site associated DNA (RAD) libraries following Baird et al. (2008) with some modifications. We used 400 ng of genomic DNA per sample. Restriction digestion was performed overnight using the restriction endonuclease High‐Fidelity (HF) SbfI (New England Biolabs). We multiplexed about 35 individuals per library after the ligation step using TruSeq P1 adapters and custom barcodes of 5–8 bp length. The libraries were sheared for 2 min on a COVARIS M220 Focused‐ultrasonicator. In the first size selection step, sheared fragments between 300 and 700 bp were size selected on a SageELF machine and in a second size selection step by using Agencourt AMPure XP beads. All libraries were single‐end sequenced on a single Illumina HiSeq 2500 lane at the Centre of Integrative Genomics (CIG), University of Lausanne. We added 12.5% bacteriophage PhiX genomic DNA to each library to increase sequence complexity.

The reads were demultiplexed and trimmed to 90 bp length using the process_radtags script from Stacks (Catchen et al., 2013), correcting single errors in the barcode and the restriction site, and discarding reads with incomplete restriction sites. Using the FastX‐toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html by Hannon Lab), all reads with a Phred quality score below 10, and reads with more than 5% of bases with a quality score below 30 were removed. The reads of each individual were aligned to the Symphodus melops reference genome of Mattingsdal et al. (2018) consisting of 5,040 contigs (N50 = 461 652 bp) and a total size of 614 Mbp. Alignment was performed end‐to‐end with Bowtie2 (Langmead & Salzberg, 2012) using default parameters. Base score recalibration was performed using PhiX reads, as described in Marques et al. (2016). Variants and genotypes were called with Genome Analysis Toolkit (GATK) Unified Genotyper version 4.0 (McKenna et al., 2010) across all species for a phylogenetic reconstruction of the Labrini tribe, and separately for all S. rostratus and S. roissali individuals for population genomic analyses. Filtering was performed with a custom script (https://github.com/joanam/scripts/raw/master/filter_vcf.py) and VCFtools version 0.1.14 (Danecek et al., 2011). Single nucleotide polymorphisms (SNPs) were required to have a minimum quality value of 30, and genotypes were required to have a minimum quality value of 30, and a minimum sequencing depth of 10 reads. Indels, sites within 5 bp of an indel and multiallelic SNPs were removed. Additionally, we used VCFtools to remove sites with more than 25% missing data, as well as sites with a mean genotype depth greater than median ×1.5 (>83×), as these are expected to be enriched for paralogs. For the VCF file containing the genomic dataset of all species used for the phylogenetic reconstruction, monomorphic sites were kept, whereas for all other analyses, only SNPs were used.

Phylogenetic tree reconstruction

The genomic dataset, including monomorphic sites, of 21 species was converted into a phylip file with a custom script (https://github.com/joanam/scripts/raw/master/vcf2phylip.py). We estimated the best substitution model using ModelFinder implemented in IQ‐TREE v2 (Minh et al., 2020). Next, we conducted a maximum‐likelihood tree inference with 1000 rapid bootstraps with the best‐fit substitution model according to BIC (K3Pu+F+R4) in IQ‐TREE. To obtain a balanced tree, we used only the five individuals with least missing data each of S. roissali and S. rostratus. The tree was subsequently mid‐point rooted using Figtree (http://tree.bio.ed.ac.uk/software/figtree/).

Behavioural data analysis

Microhabitat association

To assess if colour morphs are associated with different microhabitats, we performed observational studies of our two focal species, S. roissali and S. rostratus, in Elba and Corsica in summer and autumn of 2017, respectively (Table S1). Our data encompass observations of a total of 86 individuals of S. roissali (75 on Elba and 11 on Corsica) and 23 individuals of S. rostratus (13 on Elba and 10 on Corsica, Table S3). Sampling locations and sampling effort per species and site are shown in Figure 1.

FIGURE 1

Sampling sites and morph frequencies among fish with behavioural data. The map shows the sampling sites on Elba (site 1) and Corsica (sites 2–5) and the proportion of green, brown and intermediate individuals for which we obtained behavioural data in Symphodus rostratus (a) and Symphodus roissali (b). The number of fish observed is reflected in the size of the pie charts. Drawings by SNF

Specific mosaic microhabitat sections (rocky substrate with brown algae versus P. oceanica patches, Figure S1) were chosen to conduct observations while snorkelling. Focal fishes of S. roissali and S. rostratus were tracked irrespective of colour and observed for a maximum amount of 5 min or until visual contact was lost. The fishes were observed from a minimum distance of approximately 2 m to avoid reactions triggered by the human presence. Standard length of each fish was estimated with an underwater ruler. The time a focal fish of either species spent in Posidonia patches or over rocky and brown algae substrate was recorded with an underwater stopwatch (Casio illuminator and Casio F‐91W). Each site was assessed in the morning from 10:00 AM to 11:30 AM, at midday, from 12:00 PM to 1:30 PM and in the afternoon from 3:00 PM to 5:00 PM on subsequent days to account for possible daytime specific behaviour and changing water conditions (i.e., water clarity/visibility). Whenever possible, fishes were filmed with GoPro 3+ and GoPro hero 5 cameras and pictures were taken with GoPro 3+ and Olympus Stylus TG‐4 cameras on auto‐mode. Pictures and videos were analysed in raw format and were neither adjusted nor light or colour manipulated.

Observation data from Corsica were acquired from video footage for a maximum duration of 5 min per focal individual, filmed in late September 2017 at four different sites in Corsica (Figure 1). Standard length of individuals was estimated with an underwater ruler on site, or, if video footage was analysed (Corsica, 2017), from a reference in their environment, such as the width of a fully grown Neptune grass leaf. To avoid pseudo‐replication, we considered individual scale patterns, marks or scars to distinguish between individuals of the same colour morph. Colour morphs were assigned qualitatively as green, brown or intermediate in colouration by SNF and OS independently. Where available, pictures or video footage was used to confirm this classification (photo/video confirmation available for 77% of all S. roissali and for 87% of all S. rostratus observations) and the final classification was agreed upon by both SNF and OS. The fish were classified into green, if fins and most of the body patterning were green, whereas they were classified into brown when the patterning was mainly brown/greyish and the fins orange/yellow. Intermediates were difficult to assign due to high variability. If some of the patterning was olive/brown and the fins green, the fish was assigned as intermediate.

For the statistical analysis of effects predicting microhabitat association in both species, a linear model was calculated for each species separately, with the amount of time spent in P. oceanica (in seconds) as dependent variable, and colour, standard length, geographic location, water clarity (clear = visibility > 5 m, turbid = visibility < 5 m) and daytime as independent variables. The total amount of observation time spent in both microhabitats (time spent in P. oceanica patches plus time spent over rocky substrate) was added as an offset variable (lm(‘Time in Posidonia (s)’ ~ Colour + ‘Length(cm)’ + Site + Water_Clarity + Daytime +, data = tempdat, offset = (‘Time in Posidonia (s)’ + ‘Time in Rock (s)’)). The use of an offset variable to deal with proportion data has been successfully demonstrated in previous studies, such as Fauteux et al. (2012). To evaluate the model and relative importance of each independent variable (i.e., the contribution of each factor to the percent variation in our data explained), an analysis of variance (ANOVA) type II was performed (using the “drop1” function of the stats package in R), testing subsequently for one main effect after the other under the linear hypothesis of our model, to test which main effect is most likely to explain the overall observed variance. To test if the time spent in each microhabitat differed between conspecific morphs, we applied a nonparametric Wilcoxon rank sum test for the relative amount of time (in percentage) spent in each microhabitat to corroborate the output for significance from the linear model (wilcox.test(‘Rel_Time_in_Posidonia’[Colour==“morph1”], ‘Rel_Time_in_Posidonia’[Colour==“morph2”], paired = FALSE)). The Wilcoxon rank sum test was performed once for all observed individuals and once for individuals that could be distinguished by pattern, marks or colouration from pictures and video footage. A Tukey post hoc test (glht(linear_model, linfct = mcp(Colour=“Tukey”), test = adjusted(“Shaffer”))) was used to test for microhabitat association of each morph compared to another (green, intermediate and brown) via multiple comparisons, if the sample size of intermediates was sufficient. Additionally, for each species, we tested the microhabitat association separately for the Corsica and Elba datasets.

Interactions

For the more common species, S. roissali, we also tested for differences in behavioural interactions between the brown and green morphs in summer 2017 in Elba (Table S1). We recorded behaviours during all observed encounters with other wrasse individuals. Interactions were classified as either conspecific or heterospecific, and the nature of the behaviour was categorized into aggressive (chasing, biting, hitting) and non‐aggressive behaviour (swimming close‐by, feeding nearby) (Tables S5 and S6). For the statistical analysis, we pooled intraspecific and interspecific interactions because of the limited sample size (intraspecific and interspecific interactions show similar trends, not shown). Interaction counts of S. roissali were divided by total observation time to account for differences in observation time between individuals. To test for differences in aggression between conspecific morphs, we calculated a linear model, including aggressive interactions as a dependent variable, with colour morph as an independent variable and observation time as an offset variable (lm(‘Aggressive Interactions’ ~ Colour, data = Observations, offset = (‘Observation Time’)). We also provide qualitative descriptions for other, less common types of interactions.

Morphological differences between colour morphs

We tested for differences in standard length between colour morphs (Table S4) with a Wilcoxon rank sum test (wilcox.test(‘Length[cm]’[Colour==“morph1”], ‘Length[cm]’ Colour==“morph2”], paired = FALSE)). Whenever feasible, we visually identified the sex of the individuals collected in Corsica in spring 2018 (Table S1) to check for potential sex‐associated colour differences.

Next, we assessed body shape variation between the morphs in both species. Following Claverie and Wainwright (2014), we computed elongation (body depth divided by standard length). We set 10 digital landmarks (Figure S2) on images of fish collected but mostly released again between 2011 and spring 2018 in Calvi, Corsica (Table S1). Linear distances were measured with a digital calliper to the nearest millimetre (mm). The dataset constitutes images of a total of 44 individuals of S. rostratus (11 green, 4 intermediate and 29 brown) and 15 individuals of S. roissali (4 green, 1 intermediate and 10 brown). Intermediately coloured individuals were excluded from further analysis because they were too few to be analysed on their own. Images were converted to tps files by using tpsUtil (32‐bit version) and landmarks were placed in tpsDig2 (32‐bit version) (software by Rohlf, 2015). The measurement between the tip of the premaxillary and the hypural bones (commonly used to infer standard length) was affected by bending of the tail. Therefore, we used the beginning of the caudal peduncle (middle position between the end of the dorsal and anal fin base) instead of the hypural bones as the caudal point (Figure S2). The distance between landmarks (2) and (7) (body depth at the first spine of the dorsal fin) was used to infer body depth.

To calculate elongation, the tps file was loaded into R, scaled for size, and linear measurements were extracted between landmarks (1) and (4) for standard length and between landmarks (2) and (7) for body depth. Body depth was divided by standard length to obtain elongation. Landmarks other than 1/4 and 2/7 were used to correct for potential body displacements (e.g., mouth opening, bending). We calculated a linear model to test for differences in elongation between colour morphs (lm(Elongation ~ Colour, data = Morpho_data)).

Genomic differentiation between colour morphs of S. rostratus and S. roissali

All samples of S. roissali and S. rostratus used for RAD sequencing were collected between 2011 and 2018 at the STARESO research station in Calvi, Corsica (Tables S1 and S3). To investigate potential genome‐wide differentiation between the green and brown morphs in S. rostratus and in S. roissali, we performed a genomic principal component analysis (PCA) based on SNPs for each species separately using the R package SNPRelate (Zheng et al., 2012). We extracted SNPs for each species using a minor allele frequency cut‐off of 5% and performed linkage disequilibrium (LD) pruning with a custom script (https://github.com/joanam/scripts/ldPruning.sh) using an R2 threshold of 0.3 in 50 kb windows, shifting by 10 kb. Sites were required to be sequenced in at least five individuals per morph for S. rostratus (of eight green and 14 brown) and in at least two individuals per morph for S. roissali (of four green and seven brown).

Next, we tested for genome‐wide differentiation between colour morphs within S. rostratus (we had too few green S. roissali in our dataset to perform the same test for this species). F ST between green and brown morphs was estimated with Arlequin version 3.5 (Excoffier & Lischer, 2010) accounting for intra‐individual variation and allowing for maximum 20% missing data per site. We tested for significance with 100 permutations at a significance level of 5%. The matrix of genetic distances with F ST estimates per site was plotted using a sliding window approach with a minimum of five SNPs per window and a window step size of two SNPs.

Hybridization analyses

To test for gene flow between species of the Symphodus radiation, we tested for excess allele sharing with the f‐branch statistic (Malinsky et al., 2021). We pruned the phylogeny to one individual per species and used the clade including all Labrus species, Tautogolabrus adspersus and Ctenolabrus rupestris and Acantholabrus palloni, as the outgroup.

With this phylogeny, we computed f‐branch statistics for all Symphodus species with Dsuite v. 0.4 (Malinsky et al., 2021) using the filtered vcf file with all individuals. First, we computed the f‐statistics with Dsuite Dtrios specifying our tree to get f‐statistics with all trios matching the tree structure, i.e. all sister taxa compared against all more distant taxa. Next, we used DsuiteFbranch with the option “‐Z True” to output the Fbranch matrix of all comparisons and the z scores for each Fbranch value. Non‐significant values were set to 0 with a custom R script https://github.com/joanam/scripts/blob/master/removeNonsignDsuite.r. Lastly, we used dtools.py of the Dsuite package to plot the Fbranch values next to the tree.

To test for gene flow in the Labrus radiation, we computed D statistics with ADMIXTOOLS (Patterson et al., 2012) using different outgroups (see Table S9). As we do not have an outgroup to all wrasses in our phylogeny, we used f4 tests to test for gene flow between the Labrus and Symphodus with ADMIXTOOLS. We converted the vcf file to eigenstrat format with a custom bash script (https://github.com/joanam/scripts/blob/master/convertVCFtoEigenstrat.sh) using vcftools and convert f from ADMIXTOOLS. Next, we ran qpDstat with the option “printsd: YES” specified in the par file and for f4 tests “f4mode: YES.” To visualize the results, we pruned the species phylogeny to one individual per species with drop.tips() in the R package ape (Paradis & Schliep, 2019).

RESULTS

Phylogenetic reconstruction

The final genomic dataset used for the phylogenetic tree (based on 35,861polymorphic sites) contained 58 individuals belonging to 21 species of 24 known species in the tribe Labrini: A. palloni (1 sample), S. ocellatus (4), S. melops (2), S. roissali (5), S. rostratus (5), Symphodus tinca (5), Symphodus cinereus (3), Symphodus trutta (2), Symphodus bailloni (2), Symphodus caeruleus (1), Symphodus doderleini (1), Symphodus mediterraneus (5), Symphodus melanocercus (4), Centrolabrus exoletus (4), Tautoga onitis (1), T. adspersus (2), C. rupestris (3), Labrus mixtus (2), Lavandula viridis (2), Labrus merula (3) and L. bergylta (1). The only missing species are the Mediterranean wrasse Lapanella fasciata, the West African Lapanella guineensis, and the North American Tautogolabrus brandaonis. All species were supported by bootstrap support values of at least 90% (Figure 2). The two radiations (i.e., genera Labrus and Symphodus) are clearly separated. S. melanocercus groups with C. exoletus, as was predicted based on phenotype (Almada et al., 2003). The Macaronesian endemics caeruleus and trutta, formerly assigned to Centrolabrus, belong to the genus Symphodus and are closely related to the East Atlantic S. melops and very close to each other.

FIGURE 2

Phylogenomic tree of the Labrini tribe. Maximum‐likelihood tree of 21 Labrini species based on 669 324 sites, including 65 708 SNPs (single nucleotide polymorphisms). Bootstrap support estimates above 90 are shown with circles on each node. Species with green‐brown polymorphism are indicated with bold labels and depicted on the right. The asterisk denotes individuals with missing data that exceed 30%. The two radiations (i.e., Labrus and Symphodus genera) are clearly separated. The polymorphic species of the Symphodus radiation are not each other's closest relatives, unlike polymorphic Labrus species, which are all very closely related to each other. The recent radiations are labelled with arrows. The photos of Symphodus rostratus, Symphodus roissali, Symphodus ocellatus, Labrus merula and Labrus viridis were taken by OS. The photos of Labrus bergylta, and of the brown and the green Symphodus cinereus were taken by Paul Newland, Fernando Quaranta and Patrick Louisy

Notably, the colour‐polymorphic species all belong to two clades with very short internal branches, indicating that they represent rapid radiations. Taxa in lineages that show longer internal branches, potentially indicating slower speciation rates, do not show CP. The species without CP are all brown or brown with additional reddish or blueish colours, but none is green. In the genus Labrus, the three polymorphic species form a single monophyletic group, whereas in the Symphodus radiation, the polymorphic species are phylogenetically dispersed.

Behavioural data

For S. roissali, the linear model rendered colour as a highly significant predictor for the fraction of time spent in P. oceanica, while neither standard length, geographic location, water clarity nor daytime predicted microhabitat association. Among the three colour morphs, the difference between the green and brown (p = 4.78 × 10−07) and between intermediate and brown (p = 0.00012) best explained the variation in microhabitat association. The ANOVA type II analysis rendered colour as the main predictive factor for microhabitat association (p = 1.59 × 10−07), whereas none of the other factors was significant when tested as a main predictor (standard length: p = 0.10, geographic location: p = 0.29, water clarity: p = 0.88, daytime: p = 0.22). The multiple comparison with a Tukey post hoc test corroborated an association between microhabitat and colour (green vs. brown: p = 1.0 × 10−04 and intermediate vs. brown: p = 3.56 × 10−04), whereas green and intermediately coloured individuals did not differ in microhabitat association (p = 0.78). Similarly, the simpler Wilcoxon rank sum test for colour morph as predictor of microhabitat association revealed that green individuals and intermediates spent much more time in P. oceanica than brown individuals (p = 2.97 × 10−10 for green and p = 3.28 × 10−06 for intermediately coloured individuals (Figure 3a)).

FIGURE 3

Microhabitat associations of colour morphs. Boxplots showing time spent in Posidonia as percentage of total observation time, with 0 indicating that 100% of the time was spent over rocks with brown algae and 100 indicating that 100% of observed time was spent in Posidonia. Individual values are shown as dots. (a) Green (N = 18 individuals) and intermediately coloured individuals (N = 9) of Symphodus roissali spend more time in Posidonia compared to the brown morph (N = 59, both comparisons p ≤ 0.001). (b) The green morph (N = 8) of Symphodus rostratus spends more time in Posidonia compared to the brown morph (N = 13, p ≤ 0.001). Intermediately coloured individuals (N = 2) do not differ in microhabitat association from either the green or brown morph

In S. rostratus, the Wilcoxon rank sum test revealed a difference in microhabitat association between green and brown individuals with p = 0.0002, whereas intermediately coloured individuals neither differed from green nor from brown individuals (p = 0.56 and p = 0.80, respectively) (Figure 3b). Other statistical analyses could not be performed for this species because of the smaller sample sizes.

The photos and video footage available for 77% of S. roissali and 87% of S. rostratus observations confirmed all colour assignments in the field, except for two S. roissali individuals which were reassigned to brown instead of intermediate. The microhabitat analyses repeated only with those individuals for which photos or video footage was available corroborated the microhabitat association differences between morphs found with the full dataset (Results S1).

Inter‐ and intraspecific interactions of S. roissali

Among the observed individuals of green and brown S. roissali, we documented 11 interspecific interactions and 18 intraspecific interactions. For the purpose of our analyses, we pooled these interactions with conspecific or heterospecific individuals. Each recorded aggressive interaction stemmed from a different individual. The brown morph accounted for most interactions with 19 of 29 interactions (Tables S5 and S6). However, the proportion of interactions classified as aggressive was higher in the green morph than in the brown morph (linear model, p = 0.00011, Figure 4).

FIGURE 4

Aggressive interactions of Symphodus roissali morphs towards heterospecific and conspecific individuals. (a) The proportion of aggressive interactions among all interactions per morph (counts given above the bars). (b) Aggression calculated as the number of aggressive interactions divided by the total observation time of each morph, showing that the green morph exhibits more aggressive behaviour than the brown morph (p ≤ 0.001)

Qualitative description of colour morph‐associated behaviours

Symphodus roissali

During our field observations, we noticed morph‐specific predator avoidance strategies in S. roissali. Upon encounter of a potential threat, green and intermediately coloured individuals dove into Posidonia, rendering them more difficult to observe. The brown morph displayed an additional kind of evasive behaviour not observed in green individuals. Five times during our 14 days of fieldwork in Elba, we observed a brown individual flee over a distance of up to 15 m and then flip to the side, lying sideways in the brown algae close to the rocky substrate, displaying the characteristic brown and blotchy scale pattern and reduced fin movement (“lateral tilting”). This behaviour had also been observed by OS in several sites in Corsica over many years.

Symphodus rostratus

The green morph of S. rostratus displayed a “vertical hovering” behaviour, initially observed by OS throughout several years in various locations on Corsica, and confirmed by SNF in eight instances on Elba and Corsica in 2017 and 2018, respectively. Individuals were hovering vertically (snout pointing down) among leaves of P. oceanica with reduced fin movement. Vertical hovering was displayed by green S. rostratus without confrontation with an immediate threat and individuals were seen to switch between horizontal and vertical position over time; a behaviour that may enhance crypsis of green fish when in Posidonia, as whenever it is observed it is among intermediate or brown individuals.

Morphological data

We found female and male individuals among both colour morphs in S. rostratus, as depicted in Figure S3A, with a total of eight females (one green, one intermediate and six brown) and ten males (four green, one intermediate and five brown) among our preserved samples from Corsica (spring 2018). Due to lack of green individuals in S. roissali in our preserved collection from spring 2018, and the difficulty of sexing fish reliably from photos, we were unable to confirm both sexes in the green morph of this species. However, we were able to confirm brown individuals of both sexes (Figure S3B).

Standard length

The mean standard length of green and intermediately coloured individuals of S. roissali (green: 6.72 ± 1.27 cm, intermediate: 6.70 ± 1.51 cm and brown: 7.7 ± 1.22 cm) was smaller than that of brown individuals (Wilcoxon rank sum test brown vs. green: p = 0.0057 and brown vs. intermediate: p = 0.0493) (Figure 5a). Symphodus rostratus morphs did not differ in standard length (green: 9.5 ± 1.15 cm, intermediate: 10.5 ± 2.41 cm and brown: 10.15 ± 0.89 cm, Wilcoxon rank sum test for green vs. brown, green vs. intermediate and intermediate vs. brown yielded p values of p = 0.4017, p = 0. 5947 and p = 0.9291, respectively) (Figure 5b).

FIGURE 5

Standard length and body elongation in Symphodus roissali and Symphodus rostratus morphs. (a) The brown morph of Symphodus roissali with a mean standard length of 7.7 cm (±1.22) differs in size from the green morph with a mean standard length of 6.72 cm (±1.27). Most of the observed fishes range between 6 and 9 cm, which is considered a mature size. Green and intermediately coloured individuals each contain a single outlier of 3 cm standard length, which likely represent subadults. (b) There is no significant size difference between morphs of Symphodus rostratus (with mean standard lengths of 9.5 ± 1.15 cm, 10.5 ± 2.41 cm and 10.15 ± 0.89 cm, for green, intermediate and brown individuals, respectively). (c) We found no difference in elongation between the colour morphs of Symphodus roissali (p > 0.05). (d) Elongation separates the green and brown morphs of Symphodus rostratus with p < 0.001. Higher values correspond to a deeper body relative to body length. Thus, the brown morph of Symphodus rostratus is deeper bodied compared to the green morph

Body elongation

The linear model revealed that the green and brown morphs of S. rostratus differed in body elongation (p < 0.001). The brown morph generally had a deeper body compared to the green morph (Figure 5d, Table S7). In contrast, we found no significant difference in elongation between the morphs of S. roissali (p = 0.85) (Figure 5c, Table S8).

The genomic dataset of S. roissali consisted of 21 278 SNPs including three green and eight brown individuals. In the genomic PCA, there is no separation by colour morph (Figure S4A). For S. rostratus, the genomic dataset consisted of 11 973 SNPs with eight green, six intermediates and 15 brown individuals. The PCA does not separate the individuals by colour morph (Figure S4B). For S. rostratus, the overall weighted FST between S. rostratus colour morphs was not different from 0 (F ST = 0.0000081, p = 0.52). When computing F ST values in windows of 5 SNPs shifting by 2, we discovered a single window exceeding FST estimates of permuted datasets (highest permuted F ST = 0.23). This outlier window has an F ST estimate of 0.25 and is located on scaffold SYMME_00023611 highlighted with a red arrow in Figure S5. In S. roissali, we did not have enough well‐sequenced individuals to test for genomic differentiation between morphs.

Multiple species of the Symphodus radiation show excess allele sharing with each other, indicating ongoing or past gene flow (Figure S6). The four colour‐polymorphic Symphodus species all show excess allele sharing with at least one other colour‐polymorphic species. Symphodus roissali shows excess allele sharing with the polymorphic species S. ocellatus. Symphodus rostratus shows excess allele sharing with four species. The strongest excess allele sharing is with the colour‐polymorphic species S. cinereus. It is thus possible that introgression of haplotypes underlying the CPs may contribute to the high number of species harbouring colour morphs.

Similarly, in the Labrus radiation, we find evidence for excess allele sharing between L. Merula with L. Bergylta and with L. Mixtus (Figure S7, Table S9). In addition, we find evidence for excess allele sharing deep in the phylogeny (Figure S7, Table S9). The Symphodus radiation shows excess allele sharing with T. adspersus and with the Labrus radiation. The different radiation members show the same excess allele sharing patterns, indicating that gene flow affected the base of the radiations.

DISCUSSION

Phylogeny of Labrini suggests colour polymorphism is phylogenetically widespread

To our knowledge, our phylogenomic tree is the most robust and complete phylogenetic reconstruction of the Labrini tribe including 21 of 24 described species. High bootstrap values for all nodes confirm that there are two separate radiations in the Labrini tribe. The CP is found in multiple species across the radiation, implying that multiple species either have standing variation for alleles underlying different morphs or the genetic basis for plasticity in colouration.

The polymorphic Labrus species form a single monophyletic clade suggesting a polymorphic ancestor in that clade. In contrast, the polymorphic Symphodus species are dispersed within the Symphodus clade. It is possible that the CP is also ancestral in the Symphodus radiation and was secondarily lost in monomorphic species. Alternatively, alleles underlying the polymorphism may have evolved multiple times or may have been shared between species through introgression. Given the high number of closely related species with a similar polymorphism, parallel de novo evolution of the morphs in Symphodus seems rather unlikely. However, the polymorphism may have arisen independently in the genera Symphodus and Labrus. The parallel maintenance of the morphs in many species suggests an adaptive role of the CP (Butlin et al., 2013; Jamie & Meier, 2020).

Given that the colour polymorphism (CP) is found in multiple species in the two adaptive radiations but not in any of the genera that have not radiated, the CP could also play a role in adaptive radiation. The CP could lead to increased speciation rates if the colour morphs represent precursors of incipient species or increase population size by widening the niche available to colour‐polymorphic species (Gray & McKinnon, 2006; Hugall & Stuart‐Fox, 2012; Jamie & Meier, 2020). As the genera Symphodus and Labrus represent two young radiations (3.5 Mya and 135 000 years, respectively, Hanel et al., 2002) with sympatric species that show evidence of current or past hybridization (Figures S6 and S7, Table S9), it is also possible that the presence of closely related species in the adaptive radiations facilitates the maintenance of CP (Jamie & Meier, 2020). We find that all colour‐morphic Symphodus species show evidence of past hybridization with other polymorphic species and some Labrus species show excess allele sharing with each other (Figure S7). Therefore, it is possible that if an advantageous morph evolves in one species, the underlying allele might spread across multiple species via gene flow, or if a beneficial allele is lost in one species, it might be reintroduced via gene flow from another species. Given that we find evidence for gene flow between the Symphodus and the Labrus radiation, it is possible that they share highly similar colour morphs due to introgression. Whole genome sequencing would be required to identify the genes underlying the CP in the different species and to test if the alleles share a common descent or have introgressed between species.

Contrasting microhabitat association of colour morphs suggests background matching

In both species, colour is the best predictor for microhabitat association, whereby the green morph shows strong association with P. oceanica, and the brown morph with rocky substrate. This finding is consistent with theory predicting that in heterogeneous environments (i.e., a patchy mosaic of microhabitats), differently coloured individuals should also behaviourally choose backgrounds that best match their appearance (Michalis et al., 2017; Svensson, 2017). In S. roissali and S. rostratus, the colour polymorphism and the background‐matching behaviour may represent adaptations to reduce predation in both microhabitats (Bond, 2007; Duarte et al., 2016; Gray & McKinnon, 2006; Michalis et al., 2017).

Species‐specific morph differences point to alternative strategies to enhance crypsis

The morphs do not only differ in colour but also in morphology and behaviour. In S. roissali, we found that green and intermediately coloured individuals are smaller than brown ones. Even though an association between colour and size could in principle indicate ontogenetic colour changes (Duarte et al., 2016; Sanmartín‐Villar et al., 2016; Stevens et al., 2014; Wilson et al., 2007), it is unlikely that the green or intermediate individuals in our S. roissali dataset represent juveniles, as their lengths are considered mature sizes (https://www.fishbase.se/summary/Symphodus‐roissali.html, last accessed on 23 November 2021). The size difference between morphs may thus indicate that being small is more beneficial in the Posidonia microhabitat or being large is more beneficial in rocky substrate. In the brown morph of S. roissali, we observed the behaviour of individuals flipping to their side and displaying themselves laterally while reducing fin movement, a behaviour that, to the human eye, increases crypsis against the rocky substrate.

In S. rostratus, the green morph is more elongated compared to the brown morph and the green morph exhibits the behaviour of vertical hovering among Posidonia leaves. To the human eye, an elongate green morph displaying vertical hovering resembles a drifting Posidonia leaf and may also reduce detection by visual predators. As high body depth is associated with high predation pressures in other fish species (Domenici et al., 2007), the deeper body of the brown morph could be associated with higher pressure of gape‐limited predators over rocky substrate.

In conclusion, despite harbouring a similar CP, S. roissali and S. rostratus morphs display different behaviours and morphological differences that may represent alternative strategies to escape visual predation by enhancing crypsis in their preferred microhabitats.

Morph‐specific territorial behaviour

While green individuals of S. roissali frequently displayed aggressive behaviour to both conspecifics and heterospecifics, brown individuals showed little aggression, suggesting that the green morph is more aggressive than the brown morph. Considering that Posidonia is more patchily distributed and less abundant (i.e., lower in connectivity and availability) than rocky substrate in the depth range of S. roissali, the microhabitat may qualify as a limited resource. The higher aggression observed in the green morph is thus in line with findings linking territorial behaviour to the defence of a critical resource (Mayr & Berger, 1992) which functions as a protective cover that reduces the risk of visual predation (Johnsson et al., 2004).

Multivariate differences between morphs contrast with absence of genomic differentiation

Given the multivariate differences between morphs, one might expect the evolution of assortative mating preventing mismatched combinations of colour, behaviour and morphology. However, we find no evidence of genome‐wide differentiation between our eight green and 15 brown individuals in S. rostratus. Given that this result is based on 21,278 SNPs, we should have enough power to detect genome‐wide differentiation, unless it is very weak (F ST < 0.01, Willing et al., 2012). Also in S. roissali, we find no clustering by colour in the PCA in line with random mating. However, weak differentiation between morphs in this species could remain undetected due to insufficient sample size (three green and eight brown individuals).

It is possible that some of the colour, morphological and behavioural traits are plastic and respond to the same environmental or developmental cue, preventing mismatched combinations. Such interaction of plastic colouration and habitat choice or additional traits has been observed, for instance, in the shrimp Hippolyte obliquimanus (Duarte et al., 2016) or in the ambush bug (Phymata americana, Boyle & Start, 2020). While Arigoni et al. (2002) showed that colour can change plastically in S. roisalli and S. ocellatus within weeks to months, they found no plasticity in colouration in S. rostratus. Common garden experiments would be required to identify the extent of plasticity of the different traits associated with each colour morph.

If more than one trait is genetically determined, genetic linkage must somehow be maintained among the genes coding for the many traits that differ between the morphs. If colour and the associated morphological and behavioural traits are controlled by different genes on multiple chromosomes, random segregation of chromosomes is expected to produce disadvantageous trait combinations (e.g., green individual flipping to its side over brown rocks when facing a predator) (Funk et al., 2021; Sinervo & Svensson, 2002; Svennsson, 2017). Physical linkage of the genes underlying the different traits, such as in a chromosomal inversion or a supergene, could prevent a mismatch between genetically determined traits once the variants have accumulated in a supergene region (Thompson & Jiggins, 2014). Only a single window of five SNPs showed a significantly elevated FST estimate between colour morphs in S. rostratus, but a relatively small supergene could remain undetected with sparse RAD sequencing and our relatively low number of individuals. In S. roissali, we cannot test for genomic regions of high differentiation due to too low sample size of green individuals that we obtained. Alternatively, there could be a single pleiotropic gene that affects colour, morphology and behaviour. Pleiotropy affecting such diverse traits may seem unlikely, but a comparable scenario has been found in threespine sticklebacks, where the Eda gene affects body armour and schooling behaviour (Greenwood et al., 2016).

CONCLUSIONS

The green and brown morphs shared across multiple species in the Symphodus and Labrus radiations of Mediterranean wrasses likely represent a multi‐niche polymorphism (Figure 6) as an adaptation to the mosaic habitat of rocky substrate and Posidonia sea grass patches that is widespread in the Mediterranean Sea. These wrasses are thus an ideal system to study the role of polymorphisms in mosaic habitats and how they interact with adaptive radiation.

FIGURE 6

Summary of the multi‐niche polymorphism in Symphodus roissali and Symphodus rostratus. Green and brown morphs of both species likely occupy alternative fitness peaks, Posidonia oceanica and rocky substrate, respectively. (a) In Symphodus roissali, the small, green morph hiding in Posidonia contrasts with the larger brown morph, exhibiting lateral tilting display as anti‐predator behaviour. (b) In Symphodus rostratus, the green morph is more elongated and has been observed swimming nose‐down in Posidonia, whereas the brown morph is more deep‐bodied. Drawings by SNF

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

AUTHOR CONTRIBUTIONS

OS, JIM and SNF conceived the study with suggestions from HKK and SHA. OS, SNF and SG performed behavioural observations. HKK, SHA, SNF and SG collected samples. SNF carried out the lab work. SNF and JIM analysed the data. JIM and SNF wrote the first draft of the manuscript and OS, HKK and SHA contributed to writing.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.13999.

Supplementary Material

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S9

Click here for additional data file.