Introgressive hybridization and morphological transgression in the contact zone between two Mediterranean Solea species

Abstract Hybrid zones provide natural experiments where new combinations of genotypes and phenotypes are produced. Studying the reshuffling of genotypes and remodeling of phenotypes in these zones is of particular interest to document the building of reproductive isolation and the possible emergence of transgressive phenotypes that can be a source of evolutionary novelties. Here, we specifically investigate the morphological variation patterns associated with introgressive hybridization between two species of sole, Solea senegalensis and Solea aegyptiaca. The relationship between genetic composition at nuclear loci and individual body shape variation was studied in four populations sampled across the hybrid zone located in northern Tunisia. A strong correlation between genetic and phenotypic variation was observed among all individuals but not within populations, including the two most admixed ones. Morphological convergence between parental species was observed close to the contact zone. Nevertheless, the samples taken closest to the hybrid zone also displayed deviant segregation of genotypes and phenotypes, as well as transgressive phenotypes. In these samples, deviant body shape variation could be partly attributed to a reduced condition index, and the distorted genetic composition was most likely due to missing allelic combinations. These results were interpreted as an indication of hybrid breakdown, which likely contributes to postmating reproductive isolation between the two species.


| INTRODUCTION
Secondary contacts between closely related species, resulting from either natural processes or anthropogenic activity, often lead to the formation of hybrid zones in which populations with divergent genomes have the potential to exchange genes. These zones have been shown to act as semipermeable barriers to gene flow that selectively filter introgression by preventing genomic regions involved in reproductive isolation to be exchanged among species (Barton & Gale, 1993;Barton & Hewitt, 1985;Harrison, 1993;Harrison & Larson, 2014).
Depending on the traits considered, introgressive hybridization may lead to morphological convergence of parental populations if the underlying divergent genes behave neutrally and readily introgress upon secondary contact (Grant & Grant, 2002). A simple dilution of the parental phenotypic differences whereby hybrids display intermediate phenotypes compared to the two parental species is often This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. observed in the populations located near the center of a hybrid zone (Mayr, 1963). Indeed, additive traits determined by quantitative traits loci (QTLs) with directional effects (with each QTL having an effect in one parental population but not in the other) should appear intermediate in hybrids compared to their parental populations (Albertson & Kocher, 2005;Rieseberg & Willis, 2007). Nevertheless, transgressive hybrid phenotypes that exceed the range of parental phenotypic variation have also been reported, in particular in the context of hybrid zones (Rieseberg, Archer, & Wayne, 1999;Rieseberg, Widmer, Arntz, & Burke, 2003). This phenomenon, referred to as transgressive segregation (Bell & Travis, 2005;Rieseberg, Archer et al., 1999;Stelkens, Schmid, Selz, & Seehausen, 2009), may be attributed on the contrary to traits encoded by QTLs with antagonistic effects in each parental population (Rieseberg, Archer et al., 1999).
Alternatively, morphological differences between species can persist despite genetic introgression if they are themselves involved in reproductive isolation, or if the genes that control them are tightly linked with genomic regions involved in reproductive isolation (Harrison & Larson, 2014). In this case, a link between morphological trait variation and individual fitness is expected (Arnold & Hodges, 1995). Some incompatible genotypic combinations may be eliminated, thereby preventing hybrids to reach the full range of phenotypic variation that could be potentially generated by recombining the parental genomes.
Disentangling the role of the different mechanisms that can be involved in producing morphological variation in hybrids remains fairly poorly studied in many taxa. Yet, it remains an important question for the study of hybrid zones (Gay, Crochet, Bell, & Lenormand, 2008), as it may help reveal the underlying architecture of reproductive isolation (Rieseberg, Whitton, & Gardner, 1999) as well as the origin of evolutionary novelties (Nichols et al., 2015;Parsons, Son, & Albertson, 2011).
Here, as a first step toward understanding the evolutionary constraints limiting gene flow across a species boundary, we explore body shape variation patterns along a natural hybrid zone transect between two species of sole, Solea senegalensis and Solea aegyptiaca.
Along the North African coast, S. senegalensis ranges from Senegal to Tunisia, and S. aegyptiaca occurs from Tunisia to Egypt. Few ambiguous morphological features have been described to distinguish them (Quéro, Desoutter, & Lagardère, 1986), and thus, these two taxa can be considered as cryptic species. Hence, the description of their spatial distributions relies primarily on genetics (Ouanes, Bahri-Sfar, Ben Hassine, & Bonhomme, 2011;She, Autem, Kotulas, Pasteur, & Bonhomme, 1987). The two species come into contact in a 100-km-wide zone spreading from Bizerte in the north to Cap Bon in the south. This zone encompasses the Gulf of Tunis and the Bizerte Lagoon. Introgressive hybridization occurs and was evidenced using a few genetic markers in previous studies based on allozymes (She et al., 1987) and intron length polymorphisms (Ouanes et al., 2011).
This phenomenon predominantly occurs in the large lagoon of Bizerte that appears to be the main habitat in which hybrids are found. It is however unclear whether this zone is stable or in the process of widening, because the average hybrid index observed in 2008 (41.2%; Ouanes et al., 2011) increased by three times the value measured in the same place twenty years earlier (16.1%) by She et al. (1987).
In this study, we propose to assess morphological variation patterns inside this zone and analyze the cosegregation of body shape and genetic variation. We use deviant segregation patterns at genetic markers as a way to capture possible mechanisms of selection against deleterious allelic combinations in introgressed genotypes. We also use body shape variation in admixed populations to evaluate the consequence of gene flow on the condition of introgressed individuals.

| Sampling
We collected 88 adult individuals of S. senegalensis and S. aegyptiaca

| Genetic analysis
Whole genomic DNA (40 ng/μl) was extracted from fin clips using Qiagen DNeasy blood and tissue kit. Four EPIC loci (GH2, Am2B-2, CK6-2, and Met1) that were previously used to distinguish S. senegalensis and S. aegyptiaca outside the hybrid zone by Ouanes et al. (2011) were amplified by PCR under the same conditions. PCR products were separated by electrophoresis on a 1% acrylamide gel, and individual genotypes were subsequently scored with FMBIO II (Hitachi) using an internal size standard. Genotype scoring was performed twice by two different persons and then checked for consistency.
Discriminant correspondence analysis (DCA) based on group centroids as implemented in Genetix 4.05.2 (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 1996) was used to visualize the partitioning of genetic variation among individuals. This multivariate analysis is particularly well suited to describe the genetic composition of individuals in a gradient of admixture, as often encountered in hybrid zones. Genetic differentiation estimated by F ST was assessed for each pair of samples and tested using 10,000 permutations in Genetix 4.05.2.

| Morphometric analysis
Each individual was photographed on the eyed side with a Canon Digit Ixus 95 IS using a 35 mm f/2.8 lens with a fixed focal length.
All images were digitized in TPSDig 2 (Rohlf, 2006) using 21 homolog landmarks to perform geometric morphometrics analysis based on two-dimensional type I landmarks positioned according to clear homologous anatomical features ( Figure 2). In order to remove nonshape variation, we performed a generalized Procrustes analysis using the R package Geomorph 2.1.7 (Adams & Otarola-Castillo, 2013). This transformation consists in standardizing measures to control for differences in individual body size, by translating and rotating the configuration of landmarks to minimize the sum of squared distances between homologue landmarks (Zelditch, 2004). In order to test for remaining allometric effects following Procrustes transformation, we performed a multivariate regression of individual distance to the consensus shape against size using the procD.lm() function in

Geomorph.
In order to focus on morphological variation between species, we then conducted a canonical discriminant analysis (CDA) to capture the multidimensional variation associated with body shape (Albrecht, 1980;Klingenberg, Barluenga, & Meyer, 2003;Zelditch, 2004). This method assesses the total amount of variation in body shape among groups of samples, expressed in a n-dimensions space where n is the number of groups minus one. Transformation grids produced with the thin plate spline technique were used to visualize body shape changes among sampled populations. These analyses were performed with the R packages Geomorph (Adams & Otarola-Castillo, 2013) and Morpho (Schlager, 2015).

| Combining morphometric and genetic analysis
In order to evaluate the degree of similarity between genetic and shape variation among individuals, we used the procuste() function implemented in the ade4 R package (Dray & Dufour, 2007). This method performs a simple Procrustes rotation of the hyperspaces resulting from the separate multivariate analyses to minimize the differences between the two clouds of homologous points and project them jointly in a new coordinate system (Digby & Kempton, 1987;Dray, Chessel, & Thioulouse, 2003). Finally, the significance of the concordance between genetic and morphometric data after Procrustes rotation (m²) was tested with the Procrustean randomization test (Protest) (Jackson, 1995;Peres-Neto & Jackson, 2001).

| Genetic analyses
As expected from previous studies, discriminant correspondence analysis (DCA) based on genetic data clearly separated the two samples from the periphery of the contact zone on the first axis, which explained most (50.3%) of the variation contained in the dataset.

| Geometric morphometrics
Multivariate regression between Procrustes-transformed shape coordinates and size showed no significant association (p = .12), in- Hence, along this axis, the inner samples exhibited a morphological variance exceeding that of peripheral (and hence less introgressed) samples, which constitute a transgressive pattern. On the other hand, the third axis principally opposed the Gulf of Tunis and Bizerte samples that were the furthest apart along this axis (Figure 4b). This could be considered as another kind of transgressive pattern whereby the inner samples' morphological variance also exceed that of peripheral ones, but in opposite directions.
The deformation grid obtained from the Procrustes superimposition showed that the principal body shape differences between species were associated with body height and were more specifically

| Correlation between morphometry and genetics
The Procrustes rotation of the bivariate configurations (i.e., axes 1 and 2) obtained from genetic and morphometric data analysis showed a high correlation between shape and genetic variation (m 2 = 0.67, p < .001) (Figure 5b). In the new rotated plane, the main axis of morphogenetic differentiation between species was defined by a line con-

| DISCUSSION
Our analysis of genetic variation across the hybrid zone between S. senegalensis and S. aegyptiaca is in good agreement with previous Therefore, our data support that at least part of the geographic variation in body shape is explained by the neutral introgression of additive morphological traits along a gradient of admixture. Nevertheless, we found no significant correlation between morphological and genetic variation within samples, including those from the inner part of the hybrid zone. This was expected as most variation occurs between parental populations. Moreover, we only used a few genetic markers that are extremely unlikely to be linked with the loci controlling shape variation.
Our analysis reveals a different signal of variation along the second axis of the morphological space. On this axis, the two inner populations are clearly similar to each other while being differentiated from the parental outer samples which are themselves not differentiated along axis 2. This second pattern deviates from the previously mentioned mechanism of morphological convergence as it cannot be explained by additivity of the morphological QTLs alone. Thus, it suggests the existence of more complex epistatic and pleiotropic effects due to hybridization between divergent genomes. This could be due for instance to a mechanism of hybrid breakdown acting on different types of hybrid pedigree, as the population of Bizerte Lagoon is a S. senegalensis population introgressed by S. aegyptiaca alleles, whereas it is the opposite for the Gulf of Tunis. Phenotypic effects affecting all admixed genotypes in the same direction and contrasting with the absence of such effects in both parental populations could concern some fitness-related trait reflecting the general performance of individuals, such as the condition factor. Interestingly, the deformation grid on the positive part of axis 2 shows a less convex shape of the ventral part in the inner samples which could reflect a lower condition of hybrid individuals compared to those from T A B L E 1 Genetic differentiation estimated by F ST assessed for each pair of samples Therefore, the morphogenetic correlation detected along axis 2 may reflect the fact that only hybrids combine at the same time admixed genotypes and lower condition.
Finally, a more classical pattern consistent with phenotypic transgression due to QTLs with antagonistic effects is evidenced along axis 3 of the morphospace. This effect is expected to occur in opposite directions on both sides of a hybrid zone as exemplified in Table 2. In one F I G U R E 5 Plot of the genetic and morphometric data after Procrustes rotation (a) ( : genetic data, : morphometric data; Kerkennah Islands, Gulf Of Tunis, Bizerte Lagoon, and Tabarka). Correlation (m²) between the genetic and morphometric data after Procrustes rotation (b)  Introgressed Sp2