# When European meets African honeybees (Apis mellifera L.) in the tropics: Morphological changes related to genetics in Mauritius Island (South-West Indian Ocean)

Julien Galataud, Hélène Delatte, Maéva Angélique Techer, Christophe Simiand, Preeaduth Sookar, Bernard Reynaud, Johanna Clémencet

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Julien Galataud, Hélène Delatte, Maéva Angélique Techer, Christophe Simiand, Preeaduth Sookar, et al.. When European meets African honeybees (Apis mellifera L.) in the tropics: Morphological changes related to genetics in Mauritius Island (South-West Indian Ocean). PLoS ONE, 2020, 15 (11), pp.e0242053. 10.1371/journal.pone.0242053 . hal-03126702

HAL Id: hal-03126702

https://hal.univ-reunion.fr/hal-03126702

Submitted on 1 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# When European meets African honeybees (*Apis mellifera* L.) in the tropics: Morphological changes related to genetics in Mauritius Island (South-West Indian Ocean)

Julien Galataud¹*, Hélène Delatte², Maéva Angélique Techer³, Christophe Simiand², Preeaduth Sookar⁴, Bernard Reynaud⁵, Johanna Clémencet¹*

1 Université de La Réunion, UMR PVBMT, La Réunion, France, 2 CIRAD, UMR PVBMT, 7 chemin de l'Irat, Ligne Paradis, Saint Pierre, La Réunion, France, 3 Okinawa Institute of Science and Technology, Okinawa, Japan, 4 Ministry of Agro Industry and Food Security, Agricultural Services, Reduit, Mauritius, 5 Université de La Réunion, UMR PVBMT, Saint Pierre, La Réunion, France

* julien.galataud@univ-reunion.fr (JG); johanna.clemencet@univ-reunion.fr (JC)

<figure><img src="image_3.png" /></figure>

## Abstract

The previous genetic characterization of the honeybee population of Mauritius Island (Indian Ocean) revealed an ongoing process of hybridization between the first established African subspecies *Apis mellifera unicolor* and recently imported European subspecies (*A. m. ligustica*, *A. m. carnica* and *A. m. mellifera). This context offers the rare opportunity to explore the influence of hybridization between African and European honeybees on phenotypic traits out of the case largely studied of the Africanized honeybee (hybrid between *A. m. scutellata* from South Africa and European subspecies). We thus conducted geometric morphometric analyses on forewings of 283 workers genetically characterized at 14 microsatellite loci to evaluate (1) if the morphological variability coincides well with the neutral genetic variability, (2) if hybrids exhibited rather parental, intermediate or transgressive traits, and (3) to test if fluctuating asymmetry (FA) of size and shape, as a measure of developmental stability, was elevated in hybrids (due to genetic stress) and/or European bees (due to unsuitable environment) compared to African bees. A strong concordance was found between morphological variability and neutral genetic variability, especially for wing shape, based on partial least-square analyses (PLS). However, on average, the morphology of hybrids was more similar to the African bees, potentially reflecting the dynamics and direction of introgression. Significant FA for wing size as well as wing shape was detected, suggesting the overall presence of stress during the development of the studied individuals. In contrast, the asymmetry levels do not differ according to the ancestry (African, European or hybrid) of the individuals. Therefore, if ongoing hybridization contributed to increasing the genetic and phenotypic diversity of the populations and influences its adaptive potential, developmental stressors could not be identified and their evolutionary consequences remain uncertain.

<figure><img src="image_4.png" /></figure>

### OPEN ACCESS

**Citation:** Galataud J, Delatte H, Techer MA, Simiand C, Sookar P, Reynaud B, et al. (2020) When European meets African honeybees (*Apis mellifera* L.) in the tropics: Morphological changes related to genetics in Mauritius Island (South-West Indian Ocean). PLoS ONE 15(11): e0242053.
https://doi.org/10.1371/journal.pone.0242053

**Editor:** Arnar Palsson, University of Iceland, ICELAND

**Received:** October 12, 2019

**Accepted:** October 27, 2020

**Published:** November 19, 2020

**Copyright:** © 2020 Galataud et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the manuscript and its Supporting information files.

**Funding:** This study was funded by the CIRAD, the 'Conseil Régional de La Réunion' and the European Agricultural Fund for Rural Development (EAFRD). J Galataud is a recipient of a PhD fellowship from the doctoral school of La Réunion University. The authors acknowledge the Plant Protection Platform
(3P, IBISA) where genetic experiments were conducted.

**Competing interests:** The authors have declared that no competing interests exist.

**Introduction**

Hybridization, defined as the contact between previously isolated lineages that result in viable and—to a varying degree—fertile offspring [1], has been extensively facilitated by human movements and activities, which intensify and ignore natural boundaries [2, 3]. The recent establishment of gene flows through the dissemination—voluntary or not—of individuals out of their native range may result in increased phenotypic variation, especially in closely related populations with weak or non-existent reproductive barriers. Increased phenotypic variation is typically observed among individuals as a consequence of augmented genetic variability. The import of new alleles and the consequent novel genomic combinations formed by admixture and repeated backcrossings can generate a variety of phenotypes in the hybrid offspring. These can range from parental types to intermediate or transgressive (extreme) forms, depending on both gene interactions (e.g. additive effects, over-dominance, epistasis), environmental influences on gene expression (external factors, e.g. temperatures, or internal factors, e.g. metabolism), and the relative contribution of natural selection and stochastic processes (disspersal and drift). Inbreeding and exposure to unfamiliar environmental conditions may also affect the ability of individuals to undergo stable development, through disruption of coadapted gene complexes in hybrid offspring, and as a consequence of radical environmental changes for imported genomes [4, 5]. Development instability (DI) refers to slight and random deviations to perfect left-right symmetry in structures, and is classically measured as fluctuating asymmetry (FA). In different taxa, excess of FA could be associated with (human-mediated) genetic [6–10] and environmental [11–14] stresses, and reduced adaptive value or fitness of individuals [15–17]. Thus, relating both components of morphological variation, i.e. between and within individuals variation (via FA), with neutral genetic markers permits hybrid characterization, and may provide valuable insight on the extent of hybridization as well as the mechanisms shaping hybrid populations.

The biological model we used to address these questions is the Western honeybee *Apis mellifera*, whose distribution is nowadays much influenced by beekeeping activities [18]. Four major evolutionary lineages occupying a wide natural range have been described by combining morphological and molecular approaches (A: African lineage, C and M: West and South-East European lineages, O: Middle East lineage) [19–23]. Nearly thirty *A. mellifera* subspecies and various ecotypes were identified and can supposedly all interbreed due to absent or weak reproductive isolating barriers [24–26]. The highly polyandrous mating system of the honeybee with its extensive mating flights [27], coupled with large-scale migratory beekeeping and trade-in queens, promoted gene flow and hybridization between formerly isolated honeybee populations within and outside its natural range [18]. The Africanized honeybee (AHB), a hybrid between *Apis mellifera scutellata* (native of Africa) and European-derived strains, which spread across America in a few decades and supplanted the pre-established European colonies, illustrates dramatically how human-induced hybridization can promote diversification and colonization success. The genome of the AHB is characterized by a large predominance of African ancestry [28–31], most likely due to a predisposition of African traits to South American environmental conditions and genetic and behavioral biases that benefit African alleles [32]. In addition, a positive selection for Europeans-derived alleles has also been highlighted in part of the genome [31]. The AHB therefore owes its success to the human-initiated hybridization, pre-adaptations to tropical and subtropical climate, and the adaptive exploitation of imported diversity. Among the traits that characterize the AHB, the size and shape of the wings have proven to be discriminating enough to identify it quickly and efficiently [32–35]. Moreover, wing shape of African worker bees would show greater developmental stability, which could contribute to asymmetrical gene flow and the prevalence of
African ancestry if related to better fitness (e.g., better flight skills) [36]. In short, the evolution-
ary trajectory of the American honeybee took a whole new direction following hybridization,
and the study of both genetic and phenotypic markers helps to better understanding the mech-
anisms underlying this evolutionary shift.

In this study we focused on the honeybee population of Mauritius Island, located in the South-West Indian Ocean, which also experiences ongoing hybridization between African and European lineages. It offers the opportunity to study the effects of relatively recent hybridization on traits variability in tropical and insular environmental conditions. Latest genetic studies confirmed that the initial honeybee population of Mauritius belongs to the African lineage, closely related to *A. m. unicolor* [37] subspecies, endemic of Madagascar and established in surrounding archipelagos [38–40]. Mitochondrial signatures of repeated importations of European continental subspecies in the current population have been detected, and the analyses of genetic diversity revealed the occurrence of two hybridizing gene pools, with European and African ancestries sampled in near-balanced proportions (about 55% of A lineage, 44% of European C lineage, 1% of European M lineage) [39]. In this paper, morphological variability of Mauritius honeybees was examined on forewings size and venation patterns of the worker honeybees *via* geometric morphometric methods. Geo-morphometrics are now commonly used in insect micro-taxonomy, and have been proved very effective in the identification of honeybee subspecies, ecotypes and hybrids [20, 21, 23, 41]. We then related between and among individuals variation of morphometric traits with neutral genetic structuring and the degree of hybridization of the individuals previously inferred from Bayesian analysis of micro-satellite markers [40]. The relation between hybridization on developmental instabilities, i.e. individual FA, was also examined on both forewing size and shape.

**Materials and methods**

Recent genetic work by Técher et al. [39, 40] has shown that the population of Mauritius con-
sists of two hybridizing genetic clusters: one African related to *A. m. unicolor*, whose origin
(natural or anthropogenic) is not fully elucidated, and another European from successive
imports. In order to better depict the honeybee population of Mauritius, the same bees were
used for this study. They consist of workers sampled in 2012 in apiaries with 1 to 150 colonies,
all managed by humans except for one wild colony. The bees were collected at the entrance or
inside the hives, depending on the site. All of this information is compiled in Table 1.

**Genetic data**

The use of mitochondrial markers (length and sequence of the COI-COII intergenic region [20, 42, 43]) allowed Técher et al. [39] to identify the maternal origin of each individual. Three evolutionary lineages are thus represented in the sample (Table 1): the African lineage A grouping the African subspecies, the European lineage C including subspecies naturally distributed east and south of the Alps, and the European lineage M including subspecies distributed in north and west of Europe [19]. On the other hand, a Bayesian approach involving 14 microsatellite loci (A113, A24, AP55, A88, A28, A29, AP289, AP273, (A) B124, A8, A35, AP33, AP66, and AP43) [40, 44] permitted to estimate the degree of hybridization of the individuals. This translates into a hybrid index (HI), which ranges from 0 for African ancestry to 1 for European ancestry. Intermediate values reflect hybridization to varying degrees. For convenience, HI was used as a continuous variable or transformed into a discrete variable, depending on the hypothesis tested (see Table 1 and analyzes below). The raw data containing the allele frequencies of each microsatellite marker for each individual was also retrieved for analyzes.
<table><caption>Table 1. Genetic categorization and contextual data on honeybees (Apis mellifera) sampled in Mauritius.</caption><thead><tr><th rowspan="2">Site</th><th colspan="3">Maternal lineage (mtDNA)</th><th colspan="3">Hybrid index (nDNA)</th><th rowspan="2">N</th><th rowspan="2">Date</th><th rowspan="2">N colonies/apiary</th><th rowspan="2">Type</th><th rowspan="2">Sampling</th><th colspan="4">Location</th></tr><tr><th>A</th><th>C</th><th>M</th><th>Afr.</th><th>Hybrid</th><th>Eur.</th><th>na</th><th>Latitude</th><th>Longitude</th><th>GPS X</th><th>GPS Y</th></tr></thead><tbody><tr><td>MUS01</td><td></td><td>1</td><td></td><td></td><td>1</td><td></td><td>1</td><td>14/11/12</td><td>1</td><td>M</td><td>E</td><td>S 20°23'00.1"</td><td>EO 57°25'16.4"</td><td>-20.38333</td><td>57.42111</td></tr><tr><td>MUS02</td><td>19</td><td></td><td></td><td>12</td><td>7</td><td></td><td>19</td><td>14/11/12</td><td>20</td><td>M</td><td>E</td><td>S 20°19'55.4"</td><td>EO 57°31'17.5"</td><td>-20.33194</td><td>57.52139</td></tr><tr><td>MUS03</td><td></td><td>1</td><td></td><td></td><td>1</td><td></td><td>1</td><td>13/11/12</td><td>1</td><td>M</td><td>E</td><td>S 20°11'45.5"</td><td>EO 57°35'34.6"</td><td>-20.19583</td><td>57.59278</td></tr><tr><td>MUS04</td><td>9</td><td></td><td></td><td>4</td><td>4</td><td>1</td><td>9</td><td>13/11/12</td><td>16</td><td>M</td><td>E</td><td>S 20°12'58.5"</td><td>EO 57°32'20.4"</td><td>-20.21611</td><td>57.53889</td></tr><tr><td>MUS05</td><td>1</td><td>2</td><td></td><td>1</td><td>2</td><td></td><td>3</td><td>13/11/12</td><td>3</td><td>M</td><td>E</td><td>S 20°14'00.0"</td><td>EO 57°29'24.9"</td><td>-20.23333</td><td>57.49</td></tr><tr><td>MUS06</td><td>9</td><td></td><td></td><td>5</td><td>3</td><td></td><td>9</td><td>15/11/12</td><td>10</td><td>M</td><td>E</td><td>S 20°14'33.0"</td><td>EO 57°28'38.0"</td><td>-20.242451</td><td>57.477089</td></tr><tr><td>MUS07</td><td>12</td><td>1</td><td></td><td>3</td><td>9</td><td></td><td>13</td><td>15/11/12</td><td>13</td><td>M</td><td>E</td><td>S 20°14'20.7"</td><td>EO 57°27'04.1"</td><td>-20.23889</td><td>57.45111</td></tr><tr><td>MUS08</td><td>14</td><td>1</td><td>2</td><td>5</td><td>12</td><td></td><td>17</td><td>14/11/12</td><td>22</td><td>M</td><td>I</td><td>S 20°12'06.4"</td><td>EO 57°25'12.0"</td><td>-20.20167</td><td>57.42</td></tr><tr><td>MUS09</td><td>6</td><td>16</td><td>2</td><td>3</td><td>17</td><td>2</td><td>24</td><td>15/11/12</td><td>30</td><td>M</td><td>E</td><td>S 20°10'10.7"</td><td>EO 57°27'17.2"</td><td>-20.16944</td><td>57.45472</td></tr><tr><td>MUS12</td><td>19</td><td>1</td><td></td><td>9</td><td>10</td><td></td><td>20</td><td>13/11/12</td><td>23</td><td>M</td><td>E</td><td>S 20°06'00.9"</td><td>EO 57°32'27.4"</td><td>-20.1</td><td>57.54083</td></tr><tr><td>MUS13</td><td></td><td>3</td><td></td><td></td><td>3</td><td></td><td>3</td><td>12/11/12</td><td>4</td><td>M</td><td>E</td><td>S 20°06'10.5"</td><td>EO 57°34'54.3"</td><td>-20.10278</td><td>57.58167</td></tr><tr><td>MUS14</td><td>1</td><td></td><td></td><td>1</td><td></td><td></td><td>1</td><td>15/11/12</td><td>1</td><td>W</td><td>I</td><td>S 20°06'18.2"</td><td>EO 57°38'56.7"</td><td>-20.105</td><td>57.64889</td></tr><tr><td>MUS15</td><td>4</td><td>4</td><td></td><td>3</td><td>5</td><td></td><td>8</td><td>15/11/12</td><td>8</td><td>M</td><td>E</td><td>S 20°04'19.5"</td><td>EO 57°40'26.7"</td><td>-20.07194</td><td>57.67389</td></tr><tr><td>MUS16</td><td>13</td><td></td><td></td><td>8</td><td>4</td><td></td><td>13</td><td>12/11/12</td><td>15</td><td>M</td><td>E</td><td>S 20°07'00.5"</td><td>EO 57°42'03.3"</td><td>-20.11667</td><td>57.70083</td></tr><tr><td>MUS17</td><td>13</td><td>1</td><td></td><td>3</td><td>10</td><td></td><td>14</td><td>15/11/12</td><td>15</td><td>M</td><td>E</td><td>S 20°06'38.8"</td><td>EO 57°43'08.4"</td><td>-20.11056</td><td>57.71889</td></tr><tr><td>MUS18</td><td>21</td><td>1</td><td></td><td>11</td><td>11</td><td></td><td>22</td><td>16/11/12</td><td>45</td><td>M</td><td>I</td><td>S 20°07'58.5"</td><td>EO 57°43'09.5"</td><td>-20.13278</td><td>57.71917</td></tr><tr><td>MUS19</td><td>4</td><td>9</td><td></td><td>2</td><td>10</td><td>1</td><td>13</td><td>16/11/12</td><td>19</td><td>M</td><td>E</td><td>S 20°08'17.6"</td><td>EO 57°44'41.1"</td><td>-20.13806</td><td>57.74472</td></tr><tr><td>MUS20</td><td></td><td>70</td><td></td><td></td><td>48</td><td>17</td><td>70</td><td>12/11/12</td><td>150</td><td>M</td><td>E</td><td>S 20°08'56.3"</td><td>EO 57°43'35.3"</td><td>-20.14889</td><td>57.72639</td></tr><tr><td>MUS21</td><td>1</td><td>18</td><td></td><td></td><td>11</td><td>6</td><td>19</td><td>12/11/12</td><td>21</td><td>M</td><td>E</td><td>S 20°11'40.1"</td><td>EO 57°43'38.7"</td><td>-20.19444</td><td>57.72722</td></tr><tr><td>MUS22</td><td>1</td><td></td><td></td><td>1</td><td></td><td></td><td>1</td><td>13/11/12</td><td>1</td><td>M</td><td>E</td><td>S 20°14'53.3"</td><td>EO 57°40'52.9"</td><td>-20.24806</td><td>57.68111</td></tr><tr><td>MUS23</td><td>1</td><td></td><td></td><td>1</td><td></td><td></td><td>1</td><td>13/11/12</td><td>1</td><td>M</td><td>E</td><td>S 20°17'17.2"</td><td>EO 57°41'32.2"</td><td>-20.28806</td><td>57.69222</td></tr><tr><td>MUS24</td><td>2</td><td></td><td></td><td>1</td><td>1</td><td></td><td>2</td><td>16/11/12</td><td>2</td><td>M</td><td>E</td><td>S 20°25'53.2"</td><td>EO 57°39'43.4"</td><td>-20.43139</td><td>57.66194</td></tr></tbody><tfoot><tr><td>Total</td><td>150</td><td>129</td><td>4</td><td>73</td><td>169</td><td>26</td><td>15</td><td>283</td><td></td><td></td><td></td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tfoot></table>

MUS01-24, site as designated in [40]; mtDNA, mitochondrial DNA (COI-COI intergenic region); A, African lineage; C, Eastern European lineage; M, Western European lineage; nDNA, nuclear DNA; Hybrid index (HI), Bayesian estimates of individuals' admixture (see Techer et al. [40]) subdivided into three discrete categories; Afr., "pure" Africans (HI ≤ 0.1); Hybrids (0.1 < HI < 0.9); Eur., "pure" Europeans (HI ≥ 0.9) [59]; na, data not available; N, number of individuals; M, managed colony; W, wild colony; E, entrance; I; inside.

https://doi.org/10.1371/journal.pone.0242053.t001

## Wing data acquisition

The above-mentioned genetic characterizations were carried out using DNA extracted from the legs of the bees. All remaining biological material was stored in 90% ethanol at -20 °C. Under these conditions, the wings, when not damaged by capture or manipulations, kept their venation pattern intact. Thus, only individuals whose right and left wings were available and correctly preserved were used to perform morphometric measurements. For this, the wings of each individual were cut at their basis, mounted in distilled water on a micrometric slide for scaling, and then photographed with a high-resolution camera mounted on a microscope (Leica DMC5400, Leica Microsystems). Nineteen landmarks were digitized at the vein intersections following Smith et al. [45] procedure (Fig 1A). The x and y Cartesian coordinates of each landmark were collected in ImageJ version 1.49u [46] with the Point Picker plugin (http://bigwww.epfl.ch/thevenaz/point picker/) [47]. To avoid any user effect, all measurements were made by a single operator. A total of 283 individuals already genetically characterized (lineage, allele frequencies, HI) were thus characterized morphologically. In addition, 36
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig 1.** Geometric morphometrics of the forewing venation of honeybee workers (*Apis mellifera*). (A) Location of the 19 landmarks digitized on the right and (mirrored) left forewings of *Apis mellifera* workers, following Smith et al. [45] procedure. (B) Visual representation of the centroid size as computed for 19 landmarks on a forewing of *Apis mellifera*. The open circle is the centroid (i.e. the barycenter) of the landmark configuration; the segments connecting the centroid to the landmarks represent the distances used to compute centroid size. (C) Wireframe visualization of forewing shape. The landmarks are connected with straight lines depicting the wing venation pattern.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0242053.g001

individuals randomly selected were photographed and digitized twice in separate sessions to quantify measurement errors related to image acquisition and landmark positioning.

## Forewing size and shape extraction

Overall wing size was estimated by the centroid size (CS), a measure equal to the square root of the sum of the squared distances between each landmark and the centroid, i.e. the
barycenter of the landmark configuration (Fig 1B). Wing shape variables were extracted from x and y coordinates of the 19 landmarks *via* a generalized least-squares Procrustes superimposition with the following steps [48, 49] (Fig 1C): (1) right and—previously mirrored—left wings of each individual were centered on their centroid, then (2) normalized to unit centroid-size (CS = 1.0) to remove size and shape association (i.e. isometry), and (3) iteratively rotated to minimize the sum of squared distances between each configuration and the sample mean. Thus, any difference in Procrustes coordinates of corresponding landmarks between two configurations resulted from (isometry-free) shape differences, expressed in units of Procrustes distances (dimensionless). As a consequence of Procrustes superimposition steps, 4 degrees of freedom are lost (two for centering, one for scaling, and one for rotation) resulting in 38–4 = 34 shape space dimensions.

## Amounts of among individuals’ variation and asymmetry

The first step consisted in testing the presence of directional and fluctuating asymmetry while accounting for measurement errors relating to both the position and orientation of the wings in the photos and the digitization of the landmarks at vein intersections. An ANOVA on CS for size [50, 51] and a Procrustes ANOVA on Procrustes distances for shape [49, 52] were therefore conducted on the subsample of 36 individuals whose forewings were photographed and digitized twice. Both analyses follow the same two-way mixed model design and partition the total variation into components of symmetric and asymmetric variation. Following this design, the main random effect ‘individual’ represents the variation among individuals and is computed as individuals’ left-right means (corrected for any type of asymmetry). The second main and fixed effect ‘side’ represents directional asymmetry (DA) and is computed as mean of left-right differences. Fluctuating asymmetry (FA), expressed as individual left-right differences, is represented by the interaction term ‘individual × side’ (a mixed effect). Measurement errors due to imaging (ME1) and digitizing (ME2) are quantified by the variation over replicates for each individual and side. Statistical significance of the various effects was assessed with parametric Goodall’s *F* tests [53], based on summed square Procrustes distances in the case of shape [51]. Ultimately, to determine the proportion of size or shape variation explained by each factor in the whole sample, the same analyses were conducted on the complete data but in the absence of replicates, after proving the measurement errors as negligible (see Results). All the analyses mentioned above were conducted with the software MorphoJ (https://www.R-project.org/) [54].

## Congruence between morphological and neutral genetic variation

To assess the congruence between morphometric and genetic datasets, we conducted partial least squares (PLS) analyses in MorphoJ using alternatively wings size and shape (i.e., centroid sizes and Procrustes coordinates) against microsatellite allele frequencies. The PLS method is based on a singular value decomposition of the matrix of covariances between two separate blocks of variables [55, 56]. It resulted in paired axes uncorrelated and ordered by decreasing singular values, each forming linear combinations of the variables from each dataset and accounting for as much as possible of the overall covariance. The number of axes was equal to the number of dimensions in the smaller dataset (i.e. 1 for size and 34 for shape). The overall association between morphological and genetic datasets was assessed by the Escoufier’s RV coefficient [57], a multivariate generalization of the squared Pearson correlation coefficient (r). Singular values, correlations between paired PLS axes as well as RV coefficients were statistically tested by randomly permuting the individuals (rows) in the shape block (10,000 iterations), so that the covariance structure between blocks was modified without altering the
variance within blocks [55, 56]. Additionally, to parallel dominant patterns of covariance between blocks to the major patterns of variation within each block, PLS singular axes were compared to their counterparts from principal component analyses conducted separately on shape and genetic datasets. For this purpose, angles and correlations between corresponding PLS and PCA axes were computed and statistically tested in MorphoJ [58]. Regarding genetic data, the results of multivariate analyses were faced with the Bayesian clustering method by measuring correlation of PLS and PCA scores with hybrid index values.

## Effects of hybridization on wing morphology

The hybrid index (HI) was directly exploited to specifically address how the current intraspecific admixture between two genetic pools, one prior African related to *A. m. unicolor* and one European composed of recently imported subspecies (C and M lineages), affected wings size and shape of honeybees in Mauritius. As a first step, HI was used as a grouping variable to distinct hybridized individuals from non-hybridized. Three clusters were created following Vähä and Primmer [59] recommendations: the workers were classified as 'pure' Africans if HI was smaller than 0.1 ($N_{AFR} = 74$), as 'pure' Europeans if HI was larger than 0.9 ($N_{EUR} = 27$), and as hybrids ($N_H = 177$) if HI was contained between those values [22, 60]. This results in means/medians values of 0.06/0.06 for 'pure' Africans, 0.39/0.32 for hybrids, and 0.96/0.97 for 'pure' Europeans, denoting that hybrid values are slightly shifted towards the African genetic cluster. A one-way ANOVA with post-hoc tests was conducted to compare group means based on wings size. Regarding wings shape, given the unequal distribution of individuals in different clusters and the high dimensionality of the Procrustes dataset (34 dimensions), a between-group PCA was conducted, as recommended by Mitteroecker and Bookstein [61]. This allowed (1) examining the relative arrangement of the three clusters in a reduced space of two dimensions, (2) evaluating the significance of group separation, and (3) visualizing related shape changes. Pairwise mean differences and reliability of classification rates were assessed by a leave-one-out cross-validation and permutation tests (30,000 iterations). To determine if hybrids presented phenotypic transgression, the angle formed by the vectors joining the cluster means (Hybrids—Africans—Europeans: $\tilde{HAE}$) in the shape-space was computed and tested. This analyze was achieved in RStudio software (https://www.R-project.org/) [62], using the trajectory-analysis function from the package geomorph [63, 64]. Additionally, correlations between size, bgPC1 and bgPC2 scores and HI were estimated.

## Effect of hybridization on wing asymmetry

Size FA (fluctuating asymmetry) of each individual was given by the difference between the CS of the left and the right forewing. The amount of shape FA of each individual was synthetized in the form of univariate scores based on either Procrustes distances (Procrustes FA scores) or Mahalanobis distances (Mahalanobis FA scores), which take the covariance structure into account [56, 65]. To establish the reliability of these measures (delivered in MorphoJ when performing ANOVAs), intraclass correlation coefficients (2,1) were computed in the sub-dataset with repeats [66, 67]. The association between the amounts of FA and size variability was then explored with regression analyzes and permutation tests (10,000 iterations). HI was treated as a continuous, rather than a categorical variable to better infer the relationship between the level of hybridization and the amount of developmental instability, while providing precision, flexibility and statistical power compared to ANOVAs [68, 69]. In this perspective, linear and quadratic regressions were conducted to find the optimal fit for data trends, comparing F-ratios and fit indices (AIC, AICc and BIC) for model choice. Developmental instability, measured as FA, could indeed express a variety of patterns, characterized by linearity or
monotonous change towards some specific value [70], depending of the underlying mechanisms that drive it.

# Results

## Estimating measurement errors

Repeated ANOVA measurements were performed on a subsample of 36 honeybees, photographed and digitized twice. The relative contribution of measurement errors (ME) to the total variation for centroid size and shape, due to imaging and landmark digitizing steps (estimated by the mean squares (MS)), was small compared to the components of symmetric and asymmetric variation (Table 2). In particular, fluctuating asymmetry (FA), used as an estimate of developmental instability, was found significant, with *F* ratios sufficiently high to consider ME as negligible (*F* > 18 and 64, for shape and size, respectively) [71]. As a result, it was estimated that repeated measurements were unnecessary for the whole data set in order to obtain reliable and accurate estimates of FA.

## Quantifying among individuals variation and asymmetry

The ANOVAs performed on the whole sample (Table 3) indicated that morphological variation among individuals ('individuals' effect) accounts for the largest portion of the variation in the whole sample (99.27% for size and 91.02% for shape, *P* < 0.0001). Directional asymmetry ('sides' effect) was found significant only for shape (*F* = 16.69, *P* < 0.0001), despite a tiny contribution to sample variability (0.79%). Fluctuating asymmetry also contributes to only 0.72% of size variation, compared to 8.51% for wing shape.

## Congruence between morphology and neutral genetic variation

Partial least square (PLS) analyzes were used to measure the congruence of wing size (centroid size), and then wing shape variation (Procrustes coordinates), with microsatellite allele frequencies. In both cases, we found a relatively strong overall association between morphological and genetic blocks, resulting in a *RV* coefficient of 0.29 for size (*P*<sub>perm</sub> < 0.0001) and 0.32 for shape (*P*<sub>perm</sub> < 0.0001). The PLS analysis conducted between size (one dimension) and microsatellites resulted in a single pair of axes whose corresponding PLS scores are moderately and significantly correlated (*SV* = 220.98, *r* = 0.58, *P*<sub>perm</sub> < 0.0001). Note that the size vector ordered individuals according to their wings size, so that size-PLS scores and centroid sizes are interchangeable. The PLS analysis conducted between shape and microsatellites resulted in 34 pairs of PLS axes, of which only the first was found significant, taking up 90.02% of the total covariance. Microsatellite and shape PLS1 scores were also highly correlated (*SV*1 = 0.014, *r* = 0.71, *P*<sub>perm</sub> < 0.0001). When compared to the PCAs (S1 and S2 Figs) computed separately on the same data, the major shape and genetic PLS axes appear nearly identical to the corresponding first principal components (microsatellite PC1 = 23.35%, shape PC1 = 22.22%), with pairwise correlations displaying values greater than 0.98 and highly significant (*P* < 0.0001). Overall, these results suggest that the major dimension of variation between morphological and genetic blocks is equivalent to the dominant variation within the same blocks. The associated patterns illustrated in Fig 2A and 2B reveal a (essentially) one-dimensional morphogenetic continuum polarized between African and European lineages and reflecting the ongoing hybridization process. Besides, high similarities found between the hybrid index (HI) and the dominant genetic PCA and PLS axes (*r* = 0.98, *P* < 0.0001) prove multivariate analyses as accurate as the Bayesian method.
<table><caption>Table 2. Analysis of variance of forewing size and shape in 36 honeybees photographed and digitized twice.</caption><thead><tr><th colspan="7">Size (Centroid size)</th></tr><tr><th></th><th colspan="5">ANOVA</th><th></th></tr><tr><th>Effect</th><th>SS</th><th>df</th><th>%Var.</th><th>MS</th><th>Goodall's F</th><th>P</th></tr></thead><tbody><tr><td>Individuals</td><td>12209246.53</td><td>35</td><td>96.19%</td><td>348835.62</td><td>26.50</td><td>&lt; 0.0001</td></tr><tr><td>Sides</td><td>6216.62</td><td>1</td><td>0.05%</td><td>6216.62</td><td>0.47</td><td>0.4965</td></tr><tr><td>Individuals × sides</td><td>460811.61</td><td>35</td><td>3.63%</td><td>13166.05</td><td>64.80</td><td>&lt; 0.0001</td></tr><tr><td>ME1: imaging</td><td>14626.51</td><td>72</td><td>0.12%</td><td>203.15</td><td>11.33</td><td>&lt; 0.0001</td></tr><tr><td>ME2: landmark digitizing</td><td>2581.80</td><td>144</td><td>0.02%</td><td>17.93</td><td></td><td></td></tr></tbody><thead><tr><th colspan="7">Shape (Procrustes distances)</th></tr><tr><th></th><th colspan="5">Procrustes ANOVA</th><th></th></tr><tr><th>Effect</th><th>SS ×10<sup>5</sup></th><th>df</th><th>%Var.</th><th>MS ×10<sup>5</sup></th><th>Goodall's F</th><th>P</th></tr></thead><tbody><tr><td>Individuals</td><td>13403.27</td><td>1190</td><td>88.35%</td><td>11.26</td><td>9.70</td><td>&lt; 0.0001</td></tr><tr><td>Sides</td><td>119.44</td><td>34</td><td>0.79%</td><td>3.51</td><td>3.02</td><td>&lt; 0.0001</td></tr><tr><td>Individuals × sides</td><td>1381.96</td><td>1190</td><td>9.11%</td><td>1.16</td><td>18.83</td><td>&lt; 0.0001</td></tr><tr><td>ME1: imaging</td><td>150.97</td><td>2448</td><td>1.00%</td><td>0.06</td><td>2.63</td><td>&lt; 0.0001</td></tr><tr><td>ME2: landmark digitizing</td><td>114.90</td><td>4896</td><td>0.76%</td><td>0.02</td><td></td><td></td></tr></tbody></table>

df, degrees of freedom (in Procrustes ANOVA, df is that for conventional ANOVA multiplied by the shape dimension, i.e. 34); SS, sum of squares; MS, mean squares; % Var., percentage contributions of each variance component to the total variance, computed from SS values.

'Individuals': variation among individuals (i.e. symmetric component); 'Sides': directional asymmetry (DA), 'individuals × sides' interaction: fluctuating asymmetry (FA).

$$F_{\text{Ind.}} = \text{MS}_{\text{Ind.}} / \text{MS}_{\text{Ind.} \times \text{Sides}}$$

$$F_{\text{Sides}} = \text{MS}_{\text{Sides}} / \text{MS}_{\text{Ind.} \times \text{Sides}}$$

$$F_{\text{Ind.} \times \text{Sides}} = \text{MS}_{\text{Ind.} \times \text{Sides}} / \text{MS}_{\text{ME1}}$$

$$F_{\text{ME1}} = \text{MS}_{\text{ME1}} / \text{MS}_{\text{ME2}}$$

https://doi.org/10.1371/journal.pone.0242053.t002

Effects of hybridization on wing size

The one-way ANOVA conducted on centroid sizes determined that differences exist among clusters, and both Tukey and Games-Howell (equal variances not assumed) post-hoc tests found highly significant differences for all pairwise comparisons of means ($P_{adj} < 0.0001$). The

<table><caption>Table 3. Analysis of variance for both forewing size and shape of honeybee workers (Apis mellifera).</caption><thead><tr><th colspan="7">Size (Centroid size)</th></tr><tr><th></th><th colspan="5">ANOVA</th><th></th></tr><tr><th>Effect</th><th>SS</th><th>df</th><th>%Var.</th><th>MS</th><th>Goodall's F</th><th>P</th></tr></thead><tbody><tr><td>Individuals</td><td>21922690.35</td><td>301</td><td>99.27%</td><td>72832.86</td><td>137.46</td><td>&lt; 0.0001</td></tr><tr><td>Sides</td><td>471.81</td><td>1</td><td>0.002%</td><td>471.81</td><td>0.89</td><td>0.3461</td></tr><tr><td>Individuals × sides</td><td>159485.47</td><td>301</td><td>0.72%</td><td>529.85</td><td></td><td></td></tr></tbody><thead><tr><th colspan="7">Shape (Procrustes coordinates)</th></tr><tr><th></th><th colspan="5">Procrustes ANOVA</th><th></th></tr><tr><th>Effect</th><th>SS (×10<sup>5</sup>)</th><th>df</th><th>%Var.</th><th>MS ×10<sup>5</sup></th><th>Goodall's F</th><th>P</th></tr></thead><tbody><tr><td>Individuals</td><td>27146.97</td><td>10234</td><td>91.02%</td><td>2.65</td><td>10.69</td><td>&lt; 0.0001</td></tr><tr><td>Sides</td><td>140.00</td><td>34</td><td>0.47%</td><td>4.12</td><td>16.60</td><td>&lt; 0.0001</td></tr><tr><td>Individuals × sides</td><td>2538.83</td><td>10234</td><td>8.51%</td><td>0.25</td><td></td><td></td></tr></tbody></table>

df, degrees of freedom (in Procrustes ANOVA, df is that for conventional ANOVA multiplied by the shape dimension, i.e. 34); SS, sum of squares; MS, mean squares; % Var., percentage contributions of each variance component to the total variance, computed from SS values.

'Individuals': variation among individuals (i.e. symmetric component); 'Sides': directional asymmetry (DA), 'individuals × sides' interaction: fluctuating asymmetry (FA).

$$F_{\text{Ind.}} = \text{MS}_{\text{Ind.}} / \text{MS}_{\text{Ind.} \times \text{Sides}}$$

$$F_{\text{Sides}} = \text{MS}_{\text{Sides}} / \text{MS}_{\text{Ind.} \times \text{Sides}}$$

https://doi.org/10.1371/journal.pone.0242053.t003
<figure><img src="image_8.png" /><figcaption>**Fig 2.** Patterns of covariation between forewing morphology and neutral genetic variation assessed *via* partial least-square analyses (PLS). (A) Scatter plot of the PLS of the forewing size (centroid size) and neutral genetic variation (14 microsatellite allele frequencies). The analysis resulted in one significant pair of PLS axes (10,000 permutations) accounting for 100% of covariation since size is a one-dimensional variable. (B) Scatter plot of the PLS of the forewing shape (Procrustes coordinates) and neutral genetic variation (14 microsatellite allele frequencies). The first pair of PLS axes is shown, which is also the only pair found significant by the permutation test (10,000 perm.). In the graphs, each point represents one individual. The symbol and outline color represents the maternal origin of individuals (pink circle: African lineage, yellow circle: European C lineage, black diamond: European M lineage). The</figcaption></figure>
fill color symbolizes the hybrid index (HI), i.e. the hybridization gradient from African genotypes (pink) to European genotypes (yellow). Regressions lines with 95% confidence intervals (dotted lines) are shown.

https://doi.org/10.1371/journal.pone.0242053.g002

African cluster presents the smaller wings (6453 µm, sd = 145 µm, se = 17 µm, ci = 33 µm), the greatest ones corresponding to the European cluster (6778 µm, sd = 195 µm, se = 38 µm, ci = 77 µm), and hybrids occupying an intermediate place closer to the African cluster (6560 µm, sd = 176 µm, se = 13 µm, ci = 26 µm). Finally, as expected, wing size—treated as continuous—is positively correlated with HI (r = 0.59, P < 0.0001) and coincides well with the continuum of introgression between African and European genetic pools (Fig 3A).

## Effects of hybridization on wing shape

The between-group PCA conducted on wings shape (left-right mean of Procrustes coordinates) resulted in two orthogonal axes representing the Euclidean distances among the three cluster means. The first axis (bgPC1), accounts for 90.69% of the variation between group means, against 9.31% for bgPC2. Cross-validated classification revealed poor overall accuracy with only 57.14% of individuals correctly classified, which is largely attributable to the strong overlap of the hybrids on other clusters (correct classification rate rises to 90.20% when hybrids are excluded). ‘Pure’ Africans and ‘pure’ Europeans appear more distinct, and each presents a single misclassified individual attributed to the opposite group (‘pure’ Afr.: 72.97%, ‘pure’ Eur.: 81.48%) (Table 4). The trajectory analysis failed to find an angle statistically different from 0° between the vector connecting the mean shapes of the three clusters ($\widehat{HAE} = 26.13°$, *P* = 0.29). Instead, it revealed significant differences of magnitudes (i.e. vectors size), suggesting intermediacy rather than transgressive phenotypes in hybrids shape with greater proximity toward ‘pure’ Africans (A—E: 0.020; A—H: 0.009; E—H: 0.013, *P* < 0.0001 for all distances). The associated shape changes as captured by bgPCA components are depicted in the form of wireframe graphs [72] in Fig 4. Notice that the joint analysis of size and shape (not shown) did not improved classification rates [61]. Again, strong association is observed between bgPC1 and HI (r = 0.70, *P* < 0.0001), while it is absent for bgPC2 (r = -0.01, *P* = 0.87) (Fig 3B).

## Effects of hybridization on fluctuating asymmetry

Single measure intraclass correlation coefficients (2, 1) are reported to establish the reliability of univariate measures of size FA (CSFA) and shape FA (based on Procrustes or Mahalanobis distances), used to infer levels of developmental instability. Results obtained from the subsample with repeats– 36 individuals, four measurements per wing and a single examiner—indicate that repeatability is high for CSFA (0.99, 95% IC [0.987, 0.996]) and Procrustes FA scores (0.91, 95% IC [0.86, 0.95]), but low for Mahalanobis FA scores (0.02, 95% IC [-0.10, 0.19]). Consequently, Mahalanobis FA scores were excluded from the following analyzes. Linear regression analyses failed to find significant association between size and CSFA (r = 0.07, *P* = 0.25). On the other hand, significant effect of size on Procrustes FA scores indicates that the smaller wings tend to be more asymmetrical (r = -0.12, *P* = 0.04). Regression residuals were therefore used as size-corrected FA values for shape. To assess the ability of HI to predict the amount of FA, two models were constructed: one linear and one quadratic (polynomial of degree two), both compared in terms of goodness of fit. As a result, adding a quadratic term did not improve the model fits (i.e. higher values of AIC, AICc and BIC), so that the simplest (linear) has been privileged. The slopes and coefficient of determination are not significantly
<figure><img src="image_9.png" /><figcaption>**Fig 3.** Relationships between forewing size (A) and shape (B), as captured by the between-group PCA, and the **hybrid index**. CS: average centroid sizes, i.e. mean of right and the left forewing computed for each individual, in µm; bgPC1: first axis of the between-group PCA conducted on individuals classified as “pure” Africans, “pure” European and hybrids. The symbol and outline color represents the maternal origin of individuals (Δ African lineage, ○ European C lineage, ◆ European M lineage) inferred from mitochondrial markers. Regressions lines with 95% confidence intervals (dotted lines) are shown.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0242053.g003
<table><caption>Table 4. Cross-validated classification of worker honeybees as pure Africans (A), pure Europeans (E) and hybrids (H) using forewing shape.</caption><thead><tr><th rowspan="2">Cluster</th><th colspan="3">classified as</th><th rowspan="2">Correct assignments (%)</th></tr><tr><th>“Pure” A</th><th>“Pure” E</th><th>Hybrids</th></tr></thead><tbody><tr><td>“Pure” A</td><td>54</td><td>1</td><td>20</td><td>72.00%</td></tr><tr><td>“Pure” E</td><td>1</td><td>22</td><td>4</td><td>81.48%</td></tr><tr><td>Hybrids</td><td>57</td><td>37</td><td>84</td><td>47.19%</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0242053.t004

different from zero for both CSFA (P = 0.84) and size-corrected Procrustes FA scores (P = 1.00), indicating no association between FA and HI (Fig 5).

## Discussion

This study aims to investigate the influence of ongoing hybridization in Mauritius Island, as evidenced by genetics [40], between an initial population related to the subspecies *Apis mellifera unicolor* from Madagascar and imported European subspecies, at the phenotype level. For this purpose, individual variation and fluctuating asymmetry in forewing size and shape of honeybee workers were quantified via geometric morphometric methods. These traits were then confronted with neutral genetic variation to assess the effects of hybridization with exotic subspecies on phenotypic variability and developmental instability.

## Effect of hybridization with European-related subspecies on wing size and shape

Results from multivariate analyses show that morphological variation coincides well with variation in microsatellite allele frequencies. The patterns of covariance for both wings size and shape, revealed by partial least-square analyses, are (essentially) one-dimensional and parallel the level of hybridization of the individuals (Fig 2). It results in an Afro-European continuum of admixture, which also emerges in PCAs as the dominant and only intelligible pattern of variation within shape and microsatellite datasets (except unstructured phenotypic plasticity) (S1 Fig). Overall, this highlights that demographic and stochastic processes, initiated by the recent successive imports of continental European subspecies in Mauritius, are still ongoing and prevail over possible selective processes. Oleska and Tofilski [22] previously reported such agreement between morphometrics and microsatellites, but in a context of admixture between two European subspecies, *A.m. carnica* and *A.m. mellifera*. Thus, these results reaffirm the usefulness of wing geometric morphometrics for hybrid detection and to infer patterns of introgression in honeybee populations. Yet, classification of workers as pure European, pure African and hybrids show that non-introgressed individuals are still well distinct morphologically, as suggested by wings size comparisons (Fig 3A) and shape based assignment tests (Table 4; Figs 3B and 4); besides, the analyses failed to find transgressive phenotypes (i.e., exhibiting unique morphological components) in hybrid descendant. If wings size clearly reflects ancestries of the two genetic pools, i.e. smaller wings in African-related individuals and larger wings in individuals related to European commercial subspecies [19, 36, 73–77], it can only be suspected concerning wing shape, which could have been further affected by secondary evolutionary or stochastic processes [78]. Ultimately, in the absence of breeding programs and new imports—officially banned since 2011 on the island [79, 80]—the honeybee population of Mauritius should evolve towards complete admixture. Unlike the Africanized honeybee (*A.m. scutellata* × various European subspecies), which is partly characterized by a shift of European towards African traits leading notably to smaller wings and unique shapes [35], the Mauritian
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Fig 4. Shape changes between the means of “pure” African, “pure” European and hybrid clusters along the first axis of the between-group PCA. (A) “Pure” Africans *versus* “pure” Europeans. (B) “Pure” Africans *versus* hybrids. (C) “Pure” Europeans *versus* hybrids. Shape changes are displayed as wireframe graphs symbolizing the pattern of forewing venation of *Apis mellifera*. Pink: “pure” African mean shape; yellow: “pure” European mean shape; black: hybrid mean shape.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0242053.g004
<figure><img src="image_13.png" /><figcaption>**Fig 5. Relationships between fluctuating asymmetry of forewing size (A) and shape (B) and the hybrid index.**
CSFA: difference between the centroid sizes of the right and the left forewing, in µm; Procrustes FA scores: Procrustes distance between the shapes of the right and the left forewing. The asterisk (*) designates FA values corrected for size effect. Each point represents one individual. The symbol and outline color represents the maternal origin of individuals (pink circle: African lineage, yellow circle: European C lineage, black diamond: European M lineage) inferred from mitochondrial markers. Regressions lines with 95% confidence intervals (dotted lines) are shown.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0242053.g005
honeybees exhibit on average larger forewings as well as shape changes as a result of “Europe-
anization” of its population. A comparison with the surrounding islands and archipelagos of
the SWIO where *A.m. unicolor* occurs, as well as with continental European sources, should
shed light on the particular evolutionary trajectory followed by the population of Mauritius.

Effects of hybridization with European-related subspecies on
developmental instability

Developmental instability, measured as fluctuating asymmetry of bilateral traits (i.e., left-right random variations), is closely related to both the inherent genomic architecture of individuals (heterozygosity, ploidy, intra- and intergenic interactions, genomic coadaptation) and various environmental stressors, which interact in complex—and still misunderstood—ways [81]. A common assumption states that hybridization, by instigating genetic incompatibilities or breaking coadapted gene complexes, can lead to increased developmental instability in organisms; conversely, genotypes may be less able to ensure stable development under unsuitable habitats [82, 83]. In Mauritius, where individuals from separate lineages have been brought into contact and hybridize, higher levels of developmental instability were therefore expected in hybrid descendant and/or in non-native genotypes from continental European subspecies. Significant FA has been detected in the whole sample, and affects markedly more wing shape (8.51% of variation) that wing size (<1% of variation); subtle directional asymmetry was also noticed (<1% of variation), but found significant only in shape (Table 3). These results corroborate previous studies reporting the contribution of asymmetry in wing variability in honeybees [84] and other insects [64, 85, 86]. If they suggest the overall presence of stress during the development of the studied individuals, the amounts of individual size FA and shape FA, on the other hand, do not differ according to ancestry (African, European or mixed) (Fig 5). The only trend observed and already reported by Łopuch and Tofilski [87], is that individuals with larger wings tend to present less shape asymmetry. The results, therefore, do not highlight the (potential) effects of hybridization and exposure to an unusual environment solely. Intrinsic factors specific to each subspecies, bottlenecks, or a globally stressful environment are all additional factors to consider, among others, as sources of developmental instabilities that can hide or overcome differences, difficult to detect in real conditions. To date, only two publications investigated the influence of hybridization between separate lineages on wing asymmetry in *A. mellifera*, and only one involved both European and African lineages. The first, conducted by Smith et al. [45] on two European subspecies of M and C lineages (*A.m. mellifera* and *A.m. carnica* respectively), observed no significant effect of hybridization on wing asymmetry. The second, conducted by Schneider et al. [36], in which the African subspecies *A.m. scutellata* and a European commercial strain (combination of several European subspecies) were crossbred to mimic the first stages of the emergence of the Africanized honeybee, found lower shape FA in African honeybees relative to European and hybrid workers. In theory, the similar amounts of FA observed in the African workers compared to European (mix of *A.m. ligustica*, *A.m. carnica* and *A.m. mellifera*) and hybrids in Mauritius may result from a more sensitive developmental system inherited from *A.m. unicolor* [81], or may reflect population isolation [88]. One alternative is that other stressors interact and supplant the effects of hybridization and unsuitable habitat. Indeed, Mauritius is highly subject to human pressures. Terrestrial ecosystems of the island have been largely transformed, sparing less than 2% of native habitats [89, 90], and the majority of land area is devoted to agriculture (53%), resulting in high degrees of human disturbances. Recent studies reported increased asymmetry in honeybees related to urbanization [91, 92], exposition to pesticides [93] and pollution rates [84]. Likewise, such anthropogenic stressors could impact the development of bees in Mauritius and deserve further research.
**Conclusion**

The Mauritian honeybee is experimenting intraspecific hybridization between *A.m. unicolor*-
related individuals (subspecies endemic of Madagascar) and European-related subspecies, fol-
lowing the latter's imports since the 19th century. The strong agreements found between wings
morphology, based on wings size and shape, and neutral genetic variation (14 microsatellite
markers) show that the process is still going on, and reaffirm that wings morphometry can be
used for hybrid detection in honeybees. Developmental instabilities, measured by individual
FA and observed in both wing size and shape, affect individuals in comparable levels whatever
their origin and level of hybridization. Consequently, wing asymmetry should influence
equally fitness and performance (if linked to) of the workers at the population level, and might
not contribute directly to trait displacements in the population, as it is suggested for the Afri-
canized honeybee [36]. Measuring more traits, including behavioral ones, would be necessary
and could help to clarify the evolutionary mechanisms driving the population. Human-medi-
ated hybridization has surely changed the evolutionary trajectory of the Mauritian honeybee.
In the shorter term, increased genetic and phenotypic variances could have improved its ability
to cope with growing anthropogenic threats, as illustrated by the recent and accidental intro-
duction of *Varroa destructor* in Mauritius [94], a mite implicated in elevated colony mortality
rates worldwide as a virus vector [95, 96]. Thus, results presented by Techer et al. [39, 40] and
here can beneficiate, as a reference point before *V. destructor* arrival, to the management of the
honeybee in Mauritius in order to preserve its ecosystem and economic roles. More widely, it
could provide better understanding on the impact of intraspecific hybridization and contact
with new parasites on evolutionary trajectories of populations.

**Supporting information**

S1 Fig. Patterns of variation within neutral genetic variation (A) and forewing shape (B) from principal component analyses (PCAs). Genetic and shape structures are inferred from microsatellite allele frequencies and Procrustes coordinates, respectively. Variables were centered but not scaled. Points represent individuals genotype (A) or forewing shape (B). The symbol and outline color represents the maternal origin of individuals (pink circle: African lineage, yellow circle: European C lineage, black diamond: European M lineage). The fill color symbolizes the hybrid index (HI), i.e. the hybridization gradient from African genotypes (pink) to European (yellow) genotypes, which is highly similar to the main genetic and shape (PC1) axes.
(EPS)

S2 Fig.
(EPS)

**Acknowledgments**

The authors are grateful to all beekeepers from Mauritius. The authors acknowledge the Plant Protection Platform (3P, IBISA) where all genetics experiments were conducted, and would like to thank P. Frouin for providing access to the Leica DMC5400 microscope. We thank the anonymous reviewers for their many insightful comments, which significantly contributed to improving the quality of the manuscript.

**Author Contributions**

Conceptualization: Julien Galataud, Hélène Delatte, Johanna Clémencet.
**Data curation:** Julien Galataud, Hélène Delatte, Maéva Angélique Techer, Johanna Clémencet.

**Formal analysis:** Julien Galataud.

**Funding acquisition:** Hélène Delatte, Bernard Reynaud, Johanna Clémencet.

**Investigation:** Julien Galataud, Hélène Delatte, Maéva Angélique Techer, Christophe Simiand, Preeaduth Sookar, Johanna Clémencet.

**Methodology:** Julien Galataud, Hélène Delatte, Johanna Clémencet.

**Project administration:** Bernard Reynaud, Johanna Clémencet.

**Resources:** Maéva Angélique Techer, Christophe Simiand, Preeaduth Sookar.

**Supervision:** Hélène Delatte, Bernard Reynaud, Johanna Clémencet.

**Validation:** Johanna Clémencet.

**Visualization:** Julien Galataud, Hélène Delatte, Johanna Clémencet.

**Writing – original draft:** Julien Galataud.

**Writing – review & editing:** Julien Galataud, Hélène Delatte, Maéva Angélique Techer, Bernard Reynaud, Johanna Clémencet.

## References

1. Woodruff D. Natural hybridization and hybrid zones. Systematic Biology. 1973; 22(3):213–218. https://doi.org/10.1111/2041-210x.12035

2. Mooney H, Cleland E. The evolutionary impact of invasive species. Proceedings of the National Academy of Sciences. 2001; 98(10):5446–5451. https://doi.org/10.1073/pnas.091093398 PMID: 11344292

3. Vallejo-Marín M, Hiscock S. Hybridization and hybrid speciation under global change. New Phytologist. 2016; 211(4):1170–1187. https://doi.org/10.1111/nph.14004 PMID: 27214560

4. Clarke G. Patterns of developmental stability of Chrysopa perla L. (Neuroptera: Chrysopidae) in response to environmental pollution. Environmental Entomology. 1993; 22(6):1362–1366. https://doi.org/10.1093/ee/22.6.1362

5. Møller A. Developmental Stability and Fitness: A Review. The American Naturalist. 1997; 149(5):916–932. https://doi.org/10.1086/286030 PMID: 18811255

6. Wayne R K, Modi W S, O'Brien S J. Morphological variability and asymmetry in the cheetah (*Acinonyx jubatus*), a genetically uniform species. Evolution. 1986; 40(1): 78–85. https://doi.org/10.1111/j.1558-5646.1986.tb05719.x PMID: 28564107

7. Debat V, Alibert P, David P, Paradis E, Auffray J-C. Independence between developmental stability and canalization in the skull of the house mouse. Proceedings of the Royal Society B: Biological Sciences. 2000; 267: 423–430. https://doi.org/10.1098/rspb.2000.1017 PMID: 10737397

8. Iguchi K, Watanabe K, Nishida M. Validity of fluctuating asymmetry as a gauge of genetic stress in ayu stocks. Fisheries Science. 2005; 71: 308–313. https://doi.org/10.1111/j.1444-2906.2005.00965.x

9. Rego C, Matos M, Santos M. Symmetry breaking in interspecific *Drosophila* hybrids is not due to developmental noise. Evolution. 2006; 60(4): 746–761. https://doi.org/10.1111/j.0014-3820.2006.tb01153.x PMID: 16739456

10. Eterovick P, Sloss B, Scalzo J, Scalzo J, Alford, R. Isolated frogs in a crowded world: Effects of human-caused habitat loss on frog heterozygosity and fluctuating asymmetry. Biological Conservation. 2016; 195: 52–59. https://doi.org/10.1016/j.biocon.2015.12.036

11. Marchand H, Paillat G, Montuire S, Butet A. Fluctuating asymmetry in bank vole populations (Rodentia, Arvicolinae) reflects stress caused by landscape fragmentation in the Mont-Saint-Michel Bay. Biological Journal of the Linnean Society. 2003; 80: 37–44

12. Sánchez-Chardi A, García-Pando M, López-Fuster M. Chronic exposure to environmental stressors induces fluctuating asymmetry in shrews inhabiting protected Mediterranean sites. Chemosphere. 2013; 93(6): 916–923. https://doi.org/10.1016/j.chemosphere.2013.05.056 PMID: 23800592
13. Lazić M, Miguel C, Crnobrnja-Isailović J, Kaliontzopoulou A. Effects of environmental disturbance on phenotypic variation/ an integrated assessment of canalization, developmental stability, modularity, and allometry in lizard head shape. The American Naturalist. 2015; 185(1): 4–58. https://doi.org/10.1086/679011 PMID: 25560552

14. Lezcano A, Quiroga M, Liberoff A, Van der Molen S. Marine pollution effects on the southern surf crab *Ovalipes trimaculatus* (Crustacea: Brachyura: Polybiidae) in Patagonia Argentina. Marine Pollution Bulletin. 2014; 91(2): 524–529. https://doi.org/10.1016/j.marpolbul.2014.09.038 PMID: 25293645

15. Goto S, Iwata H, Shibano S, Ohya K, Suzuki A, Ogawa H. Fruit shape variation in *Fraxinus mandshurica* var. *japonica* characterized using elliptic Fourier descriptors and the effect on flight duration. The Ecological Society of Japan. 2005; 20(6): 733–738. https://doi.org/10.1007/s11284-005-0090

16. Bakker T, Mazzi D, Kraak S. Broods of attractive three-spined stickleback males require greater paternal care. Journal of Fish Biology. 2006; 69(4): 1664–1177. https://doi.org/10.1111/j.1095-8649.2006.01194.x

17. Frey F, Bukoski M. Floral symmetry is associated with flower size and pollen production but not insect visitation rates in *Geranium robertianum* (Geraniaceae). Plant Species Biology. 2013; 29(3): 272–280. https://doi.org/10.1111/1442-1984.12021

18. De la Rúa P, Jaffé R, Dall'Olio R, Muñoz I, Serrano J. Biodiversity, conservation and current threats to European honeybees. Apidologie. 2009; 40(3):263–284. https://doi.org/10.1051/apido/2009027

19. Ruttner F. Biogeography and Taxonomy of Honeybees. Berlin, Heidelberg: Springer Berlin Heidelberg; 1988.

20. Meixner M, Pinto M, Bouga M, Kryger P, Ivanova E, Fuchs S. Standard methods for characterising subspecies and ecotypes of Apis mellifera. Journal of Apicultural Research. 2013; 52(4):1–28. https://doi.org/10.3896/ibra.1.52.4.05

21. Kandemir İ, Özkan A, Fuchs S. Reevaluation of honeybee (Apis mellifera) micro taxonomy: a geometric morphometric approach. Apidologie. 2011; 42(5):618–627. https://doi.org/10.1007/s13592-011-0063-3

22. Oleksa A, Tofilski A. Wing geometric morphometrics and microsatellite analysis provide similar discrimination of honey bee subspecies. Apidologie. 2015; 46(1):49–60. https://doi.org/10.1007/s13592-013-0212-y

23. Nawrocka A, Kandemir İ, Fuchs S, Tofilski A. Computer software for identification of honey bee subspecies and evolutionary lineages. Apidologie. 2017. https://doi.org/10.1007/s13592-017-0538-y

24. Kerr W, Bueno D. Natural crossing between Apis mellifera adansonii and Apis mellifera ligustica. Evolution. 1970; 24(1):145–148. https://doi.org/10.1111/j.1558-5646.1970.tb01747.x PMID: 28563015

25. Koeniger N, Koeniger G. Reproductive isolation among species of the genus Apis. Apidologie. 2000; 31(2):313–339. https://doi.org/10.1051/apido:2000125

26. Oleksa A, Wilde J, Tofilski A, Chybicki I. Partial reproductive isolation between European subspecies of honey bees. Apidologie. 2013; 44(5):611–619. https://doi.org/10.1007/s13592-013-0212-y

27. Soland-Reckeweg G. Genetic differentiation and hybridization in the honey bee (Apis mellifera L.) in Switzerland. Thesis, University of Bern. 2006

28. Hall HG. Parental analysis of introgressive hybridization between African and European honeybees using nuclear DNA RFLPs. Genetics. 1990; 125(3):611–621. PMID: 1974226

29. Zayed A, Whitfield C. A genome-wide signature of positive selection in ancient and recent invasive expansions of the honey bee *Apis mellifera*. Proceedings of the National Academy of Sciences. 2008; 105(9):3421–3426. https://doi.org/10.1073/pnas.0800107105 PMID: 18299560

30. Mikheyev A, Tin M, Arora J, Seeley T. Museum samples reveal rapid evolution by wild honey bees exposed to a novel parasite. Nature Communications. 2015; 6(1). https://doi.org/10.1038/ncomms8991 PMID: 26246313

31. Nelson R, Wallberg A, Simões Z, Lawson D, Webster M. Genomewide analysis of admixture and adaptation in the Africanized honeybee. Molecular Ecology. 2017; 26(14):3603–3617. https://doi.org/10.1111/mec.14122 PMID: 28378497

32. Schneider S, DeGrandi-Hoffman G, Smith D. The African Honeybee: Factors Contributing to a Successful Biological Invasion. Annual Review of Entomology. 2004; 49(1):351–376. https://doi.org/10.1146/annurev.ento.49.061802.123359 PMID: 14651468

33. Francoy T, Gonçalves L, De Jong D. Rapid morphological changes in populations of hybrids between Africanized and European honey bees. Genetics and Molecular Research. 2012; 11(3):3349–3356. https://doi.org/10.4238/2012.September.17.5 PMID: 23079828

34. Francoy T, Wittmann D, Drauschke M, Müller S, Steinhage V, Bezerra-Laure M, et al. Identification of Africanized honey bees through wing morphometrics: two fast and efficient procedures. Apidologie. 2008; 39(5):488–494. https://doi.org/10.1051/apido:2008028
35. Francoy T, Wittmann D, Steinhage V, Drauschke M, Muller S, Cunha D, et al. Morphometric and genetic changes in a population of Apis mellifera after 34 years of Africanization. Genetics and Molecular Research. 2009; 8(2):709–717. https://doi.org/10.4238/vol8-2kerr019 PMID: 19554770

36. Schneider S, Leamy L, Lewis L, DeGrandi-Hoffman G. The Influence of hybridization between African and European honeybees, *Apis mellifera*, on asymmetries in wing size and shape. Evolution. 2003; 57 (10):2350. https://doi.org/10.1111/j.0014-3820.2003.tb00247.x PMID: 14628923

37. Latreille P. Notice des espèces d'abeilles vivant en grande société, ou abeilles proprement dites, et description d'espèces nouvelles. Ann Mus Hist Nat. 1804; 5:161–178.

38. Techer M, Clémencet J, Simiand C, Portlouis G, Reynaud B, Delatte H. Genetic diversity of the honey-bee (*Apis mellifera* L.) populations in the Seychelles archipelago. Insect Conservation and Diversity. 2015; 9(1):13–26. https://doi.org/10.1111/icad.12138

39. Techer M, Clémencet J, Simiand C, Preeaduth S, Azali H, Reynaud B, et al. Large-scale mitochondrial DNA analysis of native honey bee *Apis mellifera* populations reveals a new African subgroup private to the South West Indian Ocean islands. BMC Genetics. 2017; 18(1). https://doi.org/10.1186/s12863-017-0520-8 PMID: 28577537

40. Techer M, Clémencet J, Simiand C, Turpin P, Garnery L, Reynaud B, et al. Genetic diversity and differ- entiation among insular honey bee populations in the Southwest Indian Ocean likely reflect old geo- graphical isolation and modern introductions. PloS One. 2017; 12(12):26. https://doi.org/10.1371/journal.pone.0189234 PMID: 29281653

41. Miguel I, Baylac M, Iriondo M, Manzano C, Garnery L, Estonba A. Both geometric morphometric and microsatellite data consistently support the differentiation of the *Apis mellifera* M evolutionary branch. *Apidologie*. 2011; 42(2):150–161. https://doi.org/10.1051/apido/2010048

42. Cornuet J, Garnery L, Solignac M. Putative origin and function of the intergenic region between COI and COII of Apis mellifera L. mitochondrial DNA. Genetics; 128(2):393–403. PMID: 1649072

43. Evans J, Schwarz R, Chen Y, Budge G, Cornman R, De la Rua P, et al. Standard methods for molecular research in Apis mellifera. Journal of Apicultural Research. 2013; 52(4):1–54. https://doi.org/10.3896/IBRA.1.52.4.11

44. Solignac M, Vautrin D, Loiseau A, Mougel F, Baudry E, Estoup A, et al. Five hundred and fifty microsatellite markers for the study of the honeybee (Apis mellifera L.) genome. Molecular Ecology Notes. 2003; 3(2): 307–311. https://doi.org/10.1046/j.1471-8286.2003.00436.x

45. Smith D, Crespi B, Bookstein F. Fluctuating asymmetry in the honey bee, *Apis mellifera*: effects of ploidy and hybridization. Journal of Evolutionary Biology. 1997; 10(4):551–574. https://doi.org/10.1046/j.1420-9101.1997.10040551.x

46. Schneider C, Rasband W, Eliceiri K. NIH Image to ImageJ: 25 years of image analysis. Nature Methods. 2012; 9(7):671–675. https://doi.org/10.1038/nmeth.2089 PMID: 22930834

47. Thévenaz P. Point Picker an interactive ImageJ plugin that allows storage and retrieval of a collection of landmarks. Biomedical Imaging Group, Swiss Federal Institute of Technology Lausanne. 2012. http://bigwww.epfl.ch/thevenaz/point picker/

48. Rohlf F, Slice D. Extensions of the Procrustes Method for the Optimal Superimposition of Landmarks. Systematic Zoology. 1990; 39(1):40. https://doi.org/10.2307/2992207

49. Klingenberg C, Barluenga M, Meyer A. Shape analysis of symmetric structures: quantifying variation among individuals and asymmetry. Evolution. 2002; 56(10):1909–1920. https://doi.org/10.1111/j.0014-3820.2002.tb00117.x PMID: 12449478

50. Palmer A, Strobeck C. Fluctuating Asymmetry: Measurement, Analysis, Patterns. Annual Review of Ecology and Systematics. 1986; 17(1):391–421.

51. Palmer A. Fluctuating asymmetry analyses: a primer. Developmental Instability: Its Origins and Evolu- tionary Implications. 1994;:335–364.

52. Klingenberg C, McIntyre G. Geometric morphometrics of developmental instability: analyzing patterns of fluctuating asymmetry with Procrustes methods. Evolution. 1998; 52(5):1363. https://doi.org/10.1111/j.1558-5646.1998.tb02018.x PMID: 28565401

53. Goodall C. Procrustes Methods in the Statistical Analysis of Shape. Journal of the Royal Statistical Soci- ety: Series B (Methodological). 1991; 53(2):285–321. https://doi.org/10.1111/j.2517-6161.1991.tb01825.x

54. Klingenberg C. MorphoJ: an integrated software package for geometric morphometrics. Molecular Ecol- ogy Resources. 2011; 11(2):353–357. https://doi.org/10.1111/j.1755-0998.2010.02924.x PMID: 21429143

55. Rohlf F, Corti M. Use of Two-Block Partial Least-Squares to study covariation in shape. Systematic Biol- ogy. 2000; 49(4):740–753. https://doi.org/10.1080/106351500750049806 PMID: 12116437
56. Zelditch M, Swiderski D, Sheets H. Geometric morphometrics for biologists: a primer. 2nd ed. Amsterdam: Elsevier- Academic Press; 2012.

57. Escoufier Y. Le Traitement des variables vectorielles. Biometrics. 1973; 29(4):751. https://doi.org/10.2307/2529140

58. Klingenberg C, Marugán-Lobón J. Evolutionary covariation in geometric morphometric data: analyzing integration, modularity, and allometry in a phylogenetic context. Systematic Biology. 2013; 62(4):591–610. https://doi.org/10.1093/sysbio/syt025 PMID: 23589497

59. Vähä J, Primmer C. Efficiency of model-based Bayesian methods for detecting hybrid individuals under different hybridization scenarios and with different numbers of loci. Molecular Ecology. 2006; 15(1):63–72. https://doi.org/10.1111/j.1365-294X.2005.02773.x PMID: 16367830

60. Chavez A, Saltzberg C, Kenagy G. Genetic and phenotypic variation across a hybrid zone between ecologically divergent tree squirrels (Tamiasciurus). Molecular Ecology. 2011; 20(16):3350–3366. https://doi.org/10.1111/j.1365-294X.2011.05184.x PMID: 21771139

61. Mitteroecker P, Bookstein F. Linear Discrimination, Ordination, and the Visualization of Selection Gradients in Modern Morphometrics. Evolutionary Biology. 2011; 38(1):100–114. https://doi.org/10.1007/s11692-011-9109-8

62. Team RC. A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria; 2016.

63. Collyer M, Adams D. Phenotypic trajectory analysis: comparison of shape change patterns in evolution and ecology. Hystrix. 2013; 40(1). https://doi.org/10.4404/hystrix-24.1-6298

64. Adams D, Otárola-Castillo E. Geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods in Ecology and Evolution. 2013; 4(4):393–399. https://doi.org/10.1111/2041-210x.12035

65. Klingenberg C, Monteiro L. Distances and directions in multidimensional shape spaces: implications for morphometric applications. Systematic Biology. 2005; 54(4):678–688. https://doi.org/10.1080/10635150590947258 PMID: 16126663

66. Shrout P, Fleiss J. Intraclass correlations: Uses in assessing rater reliability. Psychological Bulletin. 1979; 86(2):420–428. https://doi.org/10.1037/0033-2909.86.2.420 PMID: 18839484

67. Koo T, Li M. A Guideline of Selecting and Reporting Intraclass Correlation Coefficients for Reliability Research. Journal of Chiropractic Medicine. 2016; 15(2):155–163. https://doi.org/10.1016/j.jcm.2016.02.012 PMID: 27330520

68. Altman D, Royston P. The cost of dichotomising continuous variables. BMJ. 2006; 332(7549):1080. https://doi.org/10.1136/bmj.332.7549.1080 PMID: 16675816

69. Lazic S. Why we should use simpler models if the data allow this: relevance for ANOVA designs in experimental biology. BMC Physiology. 2008; 8(1):16. https://doi.org/10.1186/1472-6793-8-16 PMID: 18644134

70. Mikula O, Auffray J-C, Macholan M. Asymmetric size and shape variation in the Central European transect across the house mouse hybrid zone. Biological Journal of the Linnean Society. 2010; 101(1):13–27. https://doi.org/10.1111/j.1095-8312.2010.01490.x

71. Klingenberg C. Analyzing Fluctuating Asymmetry with Geometric Morphometrics: Concepts, Methods, and Applications. Symmetry. 2015; 7(2): 887. https://doi.org/10.3390/sym7020843

72. Klingenberg C. Visualizations in geometric morphometrics: how to read and how to make graphs showing shape changes. Hystrix, the Italian Journal of Mammalogy. 2013; 24(1):15–24. https://doi.org/10.4404/hystrix-24.1-7691

73. Rinderer T, Sylvester H, Brown M, Villa J, Pesante D, Collins A, et al. Field and simplified techniques for identifying Africanized and European honey bees. Apidologie. 1986; 17(1):33–48. https://doi.org/10.1051/apido:19860104

74. Buco S, Rinderer T, Sylvester H, Collins A, Lancaster V, Crewe R, et al. Morphometric differences between South American africanized and South African (Apis mellifera scutellata) honey bees. Apidologie. 1987; 18(3):217–222. https://doi.org/10.1051/apido:19870301

75. Rinderer T, Sylvester H, Buco S, Lancaster V, Herbert E, Collins A, et al. Improved simple techniques for identifying Africanized and European honey bees. Apidologie. 1987; 18(2):179–196. https://doi.org/10.1051/apido:19870208

76. Lobo J. Morphometric isozymic and mitochondrial variability of Africanized honeybees in Costa Rica. Heredity. 1995; 75(2):133–141. https://doi.org/10.1038/hdy.1995.116

77. Quezada-Euán J, Medina L. Hybridization between European and Africanized honeybees (Apis mellifera L.) in tropical Yucatan, Mexico. I. Morphometric changes in feral and managed colonies. Apidologie. 1998; 29(6):555–568. https://doi.org/10.1051/apido:19980608
78. Nolte A, Sheets H. Shape based assignment tests suggest transgressive phenotypes in natural sculpin hybrids (Teleostei, Scorpaeniformes, Cottidae). Frontiers in Zoology. 2005; 2(1):11. https://doi.org/10.1186/1742-9994-2-11 PMID: 15987531

79. Belmin R. Intérêt et faisabilité d'une Indication Géographique pour le miel de Rodrigues. Thesis, Université Henri Poincaré, Institut National Polytechnique de Lorraine. 2010.; https://hal.univ-lorraine.fr/hal-01883408

80. Jhumun G. Report on Honeybees from Mauritius for OIE conference (Regional Training Seminar on Diseases of Honey Bees for OIE National Focal Points for Animal Disease Notification). In: Ministry of Agro Industry & Food Security M, editor. OIE Ezulwini; Swaziland. Mauritius: Ministry of Agro Industry & Food Security; 2011. p. 7.

81. Klingenberg C. Phenotypic plasticity, developmental instability, and robustness: the concepts and how they are connected. Frontiers in Ecology and Evolution. 2019; 7. https://doi.org/10.3389/fevo.2019.00056

82. Graham J, Raz S, Hel-Or H, Nevo E. Fluctuating Asymmetry: Methods, Theory, and Applications. Symmetry. 2010; 2(2):466–540. https://doi.org/10.3390/sym2020466

83. Maestri R, Fornel R, Galiano D, de Freitas T. Niche Suitability Affects Development: Skull Asymmetry Increases in Less Suitable Areas. PLOS ONE. 2015; 10(4):e0122412. https://doi.org/10.1371/journal.pone.0122412 PMID: 25874364

84. Nunes L, de Araújo E, Marchini L. Fluctuating asymmetry in Apis mellifera (Hymenoptera: Apidae) as bioindicator of anthropogenic environments. Rev Biol Trop. 2015; 63(3):673–682. PMID: 26666124

85. Pélabon C, Hansen T. On the adaptive accuracy of directional asymmetry in insect wing size. Evolution. 2008; 62(11):2855–2867. https://doi.org/10.1111/j.1558-5646.2008.00495.x PMID: 18752614

86. Ludoški J, Djurakic M, Stâhls G, Milankov V. Patterns of asymmetry in wing traits of three island and one continental population of Merodon albifrons (Diptera, Syrphidae) from Greece. Evolutionary Ecology Research. 2012; 14(7):933–950.

87. Łopuch S, Tofilski A. The relationship between asymmetry, size and unusual venation in honey bees (Apis mellifera). Bulletin of Entomological Research. 2016; 106(3):304–313. https://doi.org/10.1017/S0007485315000784 PMID: 27241228

88. Băncilă R, Van Gelder I, Rotteveel E, Loman J, Arntzen J. Fluctuating asymmetry is a function of population isolation in island lizards. Journal of Zoology. 2010; 282(4):266–275. https://doi.org/10.1111/j.1469-7998.2010.00736.x

89. Florens F, Baider C, Martin G, Strasberg D. Surviving 370 years of human impact: what remains of tree diversity and structure of the lowland wet forests of oceanic island Mauritius?. Biodiversity and Conservation. 2012; 21(8):2139–2167. https://doi.org/10.1007/s10531-012-0304-4

90. Florens F. The biodiversity of Mauritius and impact resulting from anthropogenic activities including introduction of alien pest species such as the long-tailed macaque, Macaca fascicularis. Translational Science and Drug Discovery 9th International Conference of Society of Free Radical Research-Africa and 4th International Conference of International Association of Medical and Biomedical Researchers. 2015:42.

91. Banaszak-Cibicka W, Fliszkiewicz M, Langowska A, Żmihorski M. Body size and wing asymmetry in bees along an urbanization gradient. Apidologie. 2017; 49(3):297–306. https://doi.org/10.1007/s13592-017-0554-y

92. Leonard R, Wat K, McArthur C, Hochuli D. Urbanisation and wing asymmetry in the western honey bee (Apis mellifera, Linnaeus 1758) at multiple scales. PeerJ. 2018; 6:e5940. https://doi.org/10.7717/peerj.5940 PMID: 30533293

93. Ondo Zue Abaga N, Alibert P, Dousset S, Savadogo P, Savadogo M, Sedogo M. Insecticide residues in cotton soils of Burkina Faso and effects of insecticides on fluctuating asymmetry in honey bees (Apis mellifera Linnaeus). Chemosphere. 2011; 83(4):585–592. https://doi.org/10.1016/j.chemosphere.2010.12.021 PMID: 21190716

94. Esnault O, Meenowa D, Sookar P, Chauzat M, Delatte H. Spread and strain determination of Varroa destructor following its introduction to Mauritius and interactions with the bee louse Braula pretoriensis in honey bee colonies. Journal of Apicultural Research. 2018; 58(1):75–83. https://doi.org/10.1080/00218839.2018.1517987

95. Le Conte Y, Ellis M, Ritter W. Varroa mites and honey bee health: can Varroa explain part of the colony losses?. Apidologie. 2010; 41(3):353–363. https://doi.org/10.1051/apido/2010017

96. Nazzi F, Brown S, Annoscia D, Del Piccolo F, Di Prisco G, Varricchio P, et al. Synergistic Parasite-Pathogen Interactions Mediated by Host Immunity Can Drive the Collapse of Honeybee Colonies. PLoS Pathogens. 2012; 8(6):e1002735. https://doi.org/10.1371/journal.ppat.1002735 PMID: 22719246