A morphological trait involved in reproductive isolation between Drosophila sister species is sensitive to temperature

Abstract Male genitalia are usually extremely divergent between closely related species, but relatively constant within one species. Here we examine the effect of temperature on the shape of the ventral branches, a male genital structure involved in reproductive isolation, in the sister species Drosophila santomea and Drosophila yakuba. We designed a semi‐automatic measurement machine learning pipeline that can reliably identify curvatures and landmarks based on manually digitized contours of the ventral branches. With this method, we observed that temperature does not affect ventral branches in D. yakuba but that in D. santomea ventral branches tend to morph into a D. yakuba‐like shape at lower temperature. We found that male genitalia structures involved in reproductive isolation can be relatively variable within one species and can resemble the shape of closely related species’ genitalia through plasticity to temperature. Our results suggest that reproductive isolation mechanisms can be dependent on the environmental context.

| 7493 PELUFFO Et aL. and the "plasticity-first" model (Levis & Pfennig, 2016). A key feature of all these views is that the phenotypic change triggered by the plastic response, which allows the colonization of the new niches, is a phenocopy, that is, that the phenotypic change can be developmentally triggered by environmental variation or genetic variation interchangeably (Lafuente & Beldade, 2019). As we learn more about the genes mediating phenotypic plasticity (Gibert, 2017), it appears that similar phenotypic changes, either environmentally or genetically induced, can sometimes involve the same genetic loci.
For example, the same enhancer of the gene tan contributes to both phenotypic plasticity in Drosophila melanogaster (Gibert et al., 2016) and interspecific evolution between sister species Drosophila santomea and Drosophila yakuba with respect to abdomen pigmentation (Jeong et al., 2008).
Depending on the setting, plasticity can either accelerate, slow down, or have little effect on evolution and species divergence (Price et al., 2003). Speciation, the process through which lineages diverge and become reproductively isolated, involves the accumulation over time of barriers limiting interbreeding, including divergence in ecological niches, behavioral isolation, and genomic incompatibilities (Coyne & Orr, 2004). As early as 1844, anatomical differences in genitalia between closely related species were proposed to be an essential mechanism maintaining reproductive isolation, as the so-called "lock-and-key" hypothesis (Dufour, 1844;Masly, 2011).
In animals with internal fertilization, genitalia are the most rapidly evolving organs in terms of morphology (Eberhard, 1988), suggesting that a significant part of the speciation process involves anatomical divergence in genitalia. Alternatively, genital evolution can be a by-product of other evolutionary processes occurring within single lineages, independently of speciation (such as sexual selection), and lead to reproductive isolation as a by-product, when individuals attempt to hybridize with other lineages (Masly, 2011).
The lock-and-key hypothesis, even in species where it seems applicable, has been challenged by a variety of observations, including the facts that (1) genitalia in females do not differ as much as in males, (2) closely related species with conspicuous genital differences can still often produce hybrids, (3) males with laser-ablated genital organs can still copulate with no observed defect, and (4) genitalia morphology can be sensitive to temperature or nutrition (Andrade et al., 2005;Arnqvist & Thornhill, 1998;LeVasseur-Viens et al., 2015;Masly, 2011;Shapiro & Porter, 1989;Simmons, 2014 and references therein). It is thus possible that in some taxonomic groups interspecific differences in genital morphology do not contribute much to reproductive isolation.
To better comprehend the link between plasticity and speciation, careful examinations of particular cases are essential, and genital traits involved in reproductive isolation represent highly relevant model systems. How plastic are genitalia in general?
Surprisingly, few studies have examined genitalia after raising organisms in various conditions. In the water strider Aquarius remigis, the mosquito Aedes aegypti and the fly D. melanogaster, changes in larval crowding, nutrition conditions or temperature were found to affect adult body size but had little effect on the size of the external genitalia (Fairbairn, 2005;Shingleton et al., 2009;Wheeler et al., 1993). However, in two other species, the mosquito Anopheles albimanus and the fly Drosophila mediopunctata, the size and shape of the male intromittent organ were found to vary with rearing temperature (Andrade et al., 2005;Hribar, 1996). Overall, analysis of individuals sampled from the wild show that for a given arthropod or mammal species, the genitalia are usually more or less the same size whereas adult body size varies extensively (Dreyer & Shingleton, 2011;Eberhard et al., 1998 and references therein).
These observations are concordant with the "lock-and-key hypothesis," where male genitalia have to be of a particular size and shape to physically fit with the female genitalia. They are also explained by the "one-size-fits-all" hypothesis, where females appear to prefer males with genitalia of intermediate size (Eberhard et al., 1998).
In order to analyze and quantify the possible link between plasticity, reproductive isolation, and interspecific divergence, we chose to examine the effect of temperature on a male primary sexual trait likely involved in reproductive isolation between two Drosophila sister species, D. santomea and D. yakuba. These two species form an attractive system because their natural environment is relatively well characterized, they are known to hybridize, and one of their most remarkable morphological differences is a primary sexual trait that seems to be involved in a "lock-and-key" mechanism. D. santomea and D. yakuba diverged approximately 0.5-1 million years ago (Turissini & Matute, 2017). They can be crossed to generate fertile F1 females (Lachaise et al., 2000). D. santomea is endemic to the island of São Tomé, a volcanic island off the coast of Gabon (Lachaise et al., 2000), while D. yakuba is found in São Tomé and throughout sub-Saharan Africa (Lachaise et al., 1988(Lachaise et al., , 2000. In São Tomé, D. santomea lives in the mist forests at high elevations while D. yakuba is found in open habitats associated with human presence, mostly at low elevations (Llopart et al., 2005a(Llopart et al., , 2005b. Both species co-occur at mid-elevation, around 1,150 m, and hybrids have been found consistently in this hybrid zone since its discovery in 1999 (Comeault et al., 2016;Cooper et al., 2018;Lachaise et al., 2000;Llopart et al., 2005a). D. santomea being insular, it is thought that this species originated from a common ancestor with D. yakuba, which colonized the island about 0.5-1 million years ago (Cariou et al., 2001;Llopart et al., 2002;Turissini & Matute, 2017) and that the present co-occurence of D. santomea and D. yakuba in São Tomé reflects secondary colonization by D. yakuba from the African mainland, maybe during the last 500 years when Portuguese colonized the island (Cariou et al., 2001). Analysis of genomic and mitochondrial DNA sequences indicate that gene flow occurred between the D. santomea and D. yakuba more than 1,000 generations ago (Cooper et al., 2019;Turissini & Matute, 2017).
One reproductive isolating mechanism between D. yakuba and D. santomea involves a difference in ventral branches shape in the male genitalia and is the most conspicuous difference in male genitalia shape between the two species (Kamimura & Mitsumoto, 2012b;Yassin & Orgogozo, 2013; Figure 1). Ventral branches are located below the aedeagus (i.e., the insect phallus; Rice et al., 2019) and are only found in the D. yakuba complex, which comprises Drosophila teissieri, D. yakuba, and D. santomea (Yassin & Orgogozo, 2013).
In D. yakuba, spiny ventral branches insert inside female protective pouches during mating. In D. santomea, the male spines and female pouches are absent. These structures appear to play important roles during copulation. When mating with D. yakuba males, D. santomea females are wounded by the spines of the male ventral branches and they live shorter than females mating with conspecific males (Kamimura, 2012;Kamimura & Mitsumoto, 2012b;Matute & Coyne, 2010). Compared to D. teissieiri females, D. santomea females also survive less to interspecific copulation with D. mauritiana (Yassin & David, 2016). Moreover, Kamimura and Mitsumoto (2012b) reported that "copulating pairs of D. santomea males × D.
yakuba females dislodge readily when disturbed," suggesting that the spines may fasten genital coupling (Masly, 2011). We previously found that a major QTL on chromosome 3L contributes to the ventral branches shape difference between D. santomea and D. yakuba (Peluffo et al., 2015).
In São Tomé, the climate is very stable throughout the year, with only a 2.5°C-difference between the average daily temperature of the warmest month (March) and of the coldest one (July), and daily oscillations of about 5°C only (https://en.clima te-data.org/, www.world clim.org/bioclim). Based on temperature measurements at Monte Café (https://en.clima te-data.org/), we estimate that the average temperature in the hybrid zone of Bom Sucesso (1153 m) varies between 15.5°C and 18°C throughout the year. In the wild, D. santomea flies are thus likely developing mainly at temperatures around 18°C or lower.
In previous studies of ventral branch shape, flies were raised either at 21°C (Yassin & Orgogozo, 2013) or 25°C (Kamimura, 2012;Kamimura & Mitsumoto, 2012b;Peluffo et al., 2015). Here, we report that D. santomea males raised at 18°C develop spiny ventral branches comparable to those of D. yakuba raised at 25°C. This is a surprising example where organs potentially directly linked with F I G U R E 1 Landmark detection for ventral branches at 18°C and 25°C. For each individual, a picture of the ventral branches is taken (top panel). The contour is digitized by hand and smoothed (middle panel). The curvature along the contour is obtained by finite differences, which are iterated for refining; the resulting values of curvature are smoothed too (bottom panel). Smoothed curvature (vertical axis), measured in inverse micrometers, is plotted along the contour, starting from the leftmost point. The horizontal axis is the distance along the contour, called the curvilinear abscissa, measured in micrometers. Here both axes are normalized by size and represented in arbitrary units. In each plot, the left dashed vertical line is the automatic detection lower bound, the middle dashed line is the imputed global midpoint, and the right dashed line is the automatic detection upper bound (see methods). Red points represent peaks and therefore curvature maxima whereas blue points represent cavities and therefore curvature minima. Since Drosophila yakuba is not sensitive to temperature, we only show one characteristic shape. To understand where this genital structure is positioned within the male genitalia, see Figure 1 of Kamimura and Mitsumoto (2012b) reproductive isolation undergo a plastic modification similar to the difference between two sister species. To better characterize the morphological change in ventral branches shape, we developed a user-friendly method to quantify contour curvatures and automatically detect spines using machine learning. We used it to examine the plastic response of ventral branches development at 18°C and 25°C both in newly collected wild strains and in strains kept in the laboratory for many years.

| Fly rearing and imaging
Fly strains (Table 1) were kept at 22°C on standard yeast-cornmealagar medium in uncrowded conditions before the beginning of the experiments. For each strain, roughly 20 individuals were transferred from the 22°C stock to either 18°C or 25°C, kept for a minimum of two nonoverlapping adult generations. Adult males were 5 to 7 days old when frozen at −80°C for subsequent dissection. Dissection of genitalia was performed in 1X PBS at room temperature. Each genitalia was mounted on standard glass slides in DMHF (Dimethyl Hydantoin Formaldehyde, Entomopraxis) medium and kept overnight before imaging on an Olympus IX83 inverted station at 40×.

| Raw contour acquisition
All contours were digitized by the same person. Pictures were anonymized for manual contour acquisition so that the digitizer did not know the genotype. Digitization was skipped when the quality of the mounting was judged to be poor. For each picture, a custom ImageJ plugin was used to extract x, y coordinates (in pixels) of the contour. The plugin is designed to open all the pictures contained in a directory, allowing the user to manually draw a contour of the object of interest using the freehand tool of ImageJ. The raw contour is a series of points p 1 , p 2 ,…, p n in a two-dimensional space x, y where n is the number of points over which the contour passes (usually 500 < n < 1,000). The contour is open, and its endpoints are unimportant ( Figure 1). It is analyzed (and twice smoothed) as follows. Note: For each species, the most common name in the literature, the location, year of capture, and reference to origin of the strain are given. All lines are indicated in the same order as in Figure 2.

| Smoothed contour
The first layer of transformation is a rectangular smoothing filter over the raw contour to obtain the smoothed contour. At each point Here the contour smoothing parameter α, to be adjusted via learning, describes the proportion of points (relative to the total number of points forming the contour) to include in the smoothing.
This implies that the smoothed contour is 2αn points shorter (αn on each side) than the raw contour. The flatter the contour, the wider the circumscribed circle, the larger the radius r, and the smaller the curvature k. For each contour, the curvature profile is the curvature k j computed over p ′ j in p ′ 1 , p ′ n using its neighboring points (M, N,

| Raw curvature of the smoothed contour
) versus the curvilinear abscissa s j of p ′ j which is the sum of Euclidean distances from origin,

| Refined curvature
We then use this first raw curvature estimation as information to refine the curvature in a second pass. In this second measure, the refined curvature k ′ j is computed over an adaptive window of size a = 1 |k| for k < 0.1 and a = 10 otherwise: . This means that the curvature is computed over a larger distance where it is small (and curvature radius is large), which requires more smoothing, without losing the sharpness of curvature peak determination where the curvature is large.

| Smoothed curvature
To improve curvature signal to noise ratio, for each point p ′ j with coordinates x ′ j , y ′ j and refined curvature k ′ j , we compute the smoothed curvature k ′′ j as a weighted moving average with triangular weights: with w j = n, …, w i = n − |i − j| , …, w j− n = w j+ n = 0 and where is the smoothing parameter to adjust via learning. describes the proportion of points (relative to the total number of contour points) to include in the smoothing. This implies that the smoothed curvature contour is 2 n points shorter ( n on each side) than the smoothed contour.

| Landmark detection
Curvature around the start and end of the contour is noisy; it corresponds to a region of low curvature, at the beginning and end of the contour, outside of the region where we expect to find the five landmarks ( Figure 1). In addition, the contour digitization by the user, which tends to start at a precise point and to end in a long stroke, results in a slight left-right asymmetry in the curvature profile. After superimposing all smoothed curvature profiles, we choose to exclude the first and last 20% of the smoothed contour. We find that the axis of symmetry (midline) is at position 0.475 instead of 0.5 for a symmetric profile.
Landmarks are Bookstein's type 2 (local maxima of curvature) (Bookstein, 1992): maxima of the smoothed curvature for landmarks 1, 3, 5 and minima for landmarks 2 and 4. Having detected all minima and maxima, we first define landmark 3 as the maximum closest to the midline position, landmark 2 as the lowest minimum to the left of landmark 3 and landmark 1 as the maximum closest to landmark 2.
Following the same logic, we define landmark 4 as the lowest minimum to the right of landmark 3 and landmark 5 as the maximum closest to landmark 4. Having detected all five landmarks, we found that there can seldom be more than one maximum between landmark 2 and landmark 4. In such a situation, we allow resampling of landmark 3 to the highest maximum between landmarks 2 and 4. Finally, we exclude individuals that do not display all five landmarks after detection.

| Spine thrust measure
Having detected all five landmarks, we quantify form using a measure previously introduced (Peluffo et al., 2015), which is highly correlated to the Procrustes analysis principal component measure of interspecific form variation and which we called "spine thrust" (ST).
ST is a measure of how much spines are elevated above the central ridge of the ventral branches and is computed as: where Y L1 , Y L3 , and Y L5 are the Y coordinate of landmarks 1, 3, and 5, respectively. This measurement depends on the precise definition of X and Y axes. Here the X-axis is defined as the axis passing by landmarks 2 and 4 and oriented from 2 to 4, and with the Y-axis defined so that (X, Y) is an oriented orthonormal basis.

| Machine learning
Detection of maxima and minima is a simple feature detection that relies on the derivative of the smoothed curvature profile.
However, there are two parameters α, β, one for each smoothing filter (contour and curvature), which modulate the number and position of these detected maxima and minima. It is possible to explore a set of values for α and β such that the correlation be-

| Statistical analyses
All statistical analyses were performed using R version 3.4.3 (R Core Team 2016). We performed two different sets of statistical analyses to investigate how ST changes across species, year of collection and temperature. First, we fitted a standard multiple linear regression with species, year of collection, and temperature as numeric predictors using the standard R function lm(). We chose the best model based on the variance explained provided by the r 2 value. Table 2 presents the output of the lm() function using the R package jtools (v1.0.0; https://cran.r-proje ct.org/web/packa ges/jtool s/jtools.pdf) and its function export_summs() with "scale" and "transform.response" set to "TRUE" which scales and centers the response variable and reports standardized regression coefficients with their heteroskedasticity-robust standard errors. Second, we performed a regression tree analysis and performed cross-validation using recursive partitioning with the regression trees R package "rpart" version 4.1.13 (Therneau et al., 2018) and the associated function rpart() with the "ANOVA" method and obtained the approximate r 2 from a 10-fold cross-validation using the rsq.rpart() function. To confirm the importance of each factor on ST change, we also performed random forest regression analysis using the R package "RandomForest" version 4.6.14 and the randomForest() function in order. Both sets of statistical analyses investigate the role of predictors in explaining a significant part of the variance, multiple linear regression allows the use of interaction terms while regression trees are easier to interpret (James et al., 2013). In addition to these analyses, we systematically plot distribution of ST across predictors (Figure 2) showing individual values together with mean, standard errors (which directly inform about two-by-two statistical significance between groups), median, quartiles, and estimates of the 95% confidence interval of the medians, calculated as ± 1.58 × IQR √ n where IQR is the interquartile range and n the number of individuals for that IQR (Chambers et al., 1983).

| Spine thrust (ST) can be measured semiautomatically
We previously reported that the shape of ventral branches in D. santomea, D. yakuba, and their hybrids can be characterized with a set of five manually detected landmarks, which allows to calculate via simple arithmetic how much the lateral spines rise above the central ridge, as a quantitative value named "spine thrust" (ST), expressed in micrometers (Peluffo et al., 2015). The manual positioning of landmarks requires each point to be carefully positioned on the exact feature for the ST measure to be exact. It can introduce betweenuser and between-sample variability. In particular, the positioning of the three central landmarks can be equivocal and may differ between users.
To use a less biased approach and automate the process, we decided to develop a new measurement method that relies on manually digitized contours of the ventral branches, which are easier to define TA B L E 2 Results of linear model fitting. This best model shows the contribution of each explanatory variable, considered as a numerical value, and their interactions to the overall variance of ST in the full D. santomea, D. yakuba dataset (shown in Figure 2a,

| Learned α and β
We find that the same set of 30 D. yakuba and 31 D. santomea individuals is enough to identify optimal parameter values for α (contour smoothing) and β (curvature smoothing).
We find that with α = 0.025 and β = 0.055 we obtain r 2 = 0.91 ( Figure 4). Although a few other combinations of α and β yield the same r 2 (Figure 4), we choose this set because it is the one which applies the lowest degree of smoothing.

| Strong interspecific difference in ST
In total, with our semi-automated method (and after removing n = 71 individuals incorrectly dissected or mounted, 12% of total samples, with no apparent distribution bias), we phenotyped 684 individuals raised at 18°C or 25°C throughout their development, corresponding to four D. yakuba lines and seven D. santomea lines collected between 1998 and 2016 (Table 1). We checked all the automatically detected landmarks by eyes and found that 30 individuals were incorrectly digitized, with a few landmarks either missing or aberrantly positioned (see Figure 5 for a sample of such individuals), and we excluded these individuals (4% of 684) from subsequent analysis.
These aberrant landmark profiles were found in almost all the lines and at both temperatures, with no apparent distribution bias.  Table 2). Overall, and despite within-strain variability and sensitivity to temperature variation, we confirm a morphological difference of ventral branches between wild strains of D. santomea and D. yakuba using our semi-automatic method of form quantification based on ST (Figure 2a; Table 2).  (Figure 6).

F I G U R E 2
Ventral branches form is sensitive to temperature in Drosophila santomea. Isofemale lines arranged by year of collection (see Table 1 We find that the statistically significant effect of temperature on D. santomea is also statistically dependent on the year at which the strain was collected (p < 0.05, Table 2). In order to interpret our statistical analysis with multiple regression, we performed a 10-fold crossvalidated regression tree analysis on the full dataset (2 species, 11 strains, 584 individuals). The 10-fold cross-validated error rate is 0.3% and using an additive model of the shape ST ∼ species + years + temperature. We found that the variance in the dataset is first best partitioned by species and that temperature partitions the dataset best for strains collected in 2015 and 2016 (Figure 7, total variance explained as assessed by cross-validation r 2 is 0.77). To confirm those results, we also performed a random forest regression analysis with the same model as for the regression tree and found that the overall variance explained is r 2 = 0.74 and that the rank of importance of each independent variable is species > years > temperature. Altogether, our results show that in D. santomea, but not in D. yakuba, ventral branches are sensitive to temperature during development and that this effect is stronger in recently collected strains.

| The effect of temperature on spine thrust is as high as the effect of the major QTL between D. yakuba and D. santomea
To compare the effects of temperature and of interspecific genetic variation on ventral branch form, we used our previous QTL mapping dataset of ventral branch form between D. santomea and D. yakuba, which comprises 365 D. santomea backcross individuals (Peluffo et al., 2015).
In this previous study, all flies were reared at 25°C as we found that this temperature was optimal to rear both species. The five landmarks were placed manually on images of the ventral branches. A generalized Procrustes analysis was performed on a set of 365 backcross progeny individuals and a larger dataset including the backcross progeny, F1 hybrids, and parents. We found that, in both cases, the principal component PC1 explains an important part of the variance (58% in the full dataset and 41% in the backcross), that they are highly correlated (r 2 = 0.996) and that PC1 in the backcross is highly correlated to ST (spine thrust; r 2 = 0.841) and not to centroid size (r 2 = 0.038).
This QTL mapping study revealed that a 2.7Mb locus on chromosome 3L explains 30% of the mean species difference in ST, meaning that replacing one D. santomea allele at this locus with a D. yakuba allele leads to an increase in ST of about 3 μm (30% of 9 μm, Peluffo et al., 2015). Pooling all the D. santomea lines examined in the present study, we find that a change in the raising temperature from 18°C to 25°C leads to an increase in ST of about 3.4 μm (Figure 8).
We conclude that the effect of temperature is as high as the effect of genetic variation at the major interspecific genetic locus.

| A dataset-independent simultaneous quantification of shape and size
Our semi-automatic method, which relies on two simple layers of contour transformation adjusted by regression based learning, is fast and allows the measure of form variation through the simple outlining of ventral branches on 2D pictures. We note that in the future, progress in edge detection algorithms (which for now introduce too much error to measure with precision variations of the order of a few micrometers) might allow full automation from pictures to form quantification.
Having drawn contours, we could also have relied on Fourier based analyses. However, such methods require closed contours which in our case are difficult to draw since the base of the ventral branches is a complex structure which cannot be easily delimited from the cuticle of the ventral branches (Figure 1). In addition, our method is more suitable for contours in which very large and very small curvatures coexist. Furthermore, an important limitation of morphometrics analyses on landmark data (e.g., Procrustes principal component analysis) is that the PC values are dimensionless (Klingenberg, 2010) and may be difficult to relate to physical features. With our simple measure of ST obtained from the automatically detected landmarks, we are able to quantify and compare forms F I G U R E 3 Correlation between two sets of automatic measures from the same dataset. For the training dataset (31 Drosophila santomea 1563 and 30 Drosophila yakuba Oku), the same user digitized the same contours twice at one-month interval and spine thrust was automatically measured. Each point represents one individual. The y = x (black dashed line) and linear regression (full red line) are shown across studies. Importantly, because we deal with absolute geometric measurements, our method simultaneously analyzes shape and size, unlike most morphometric approaches (Claude, 2008;Klingenberg, 2010). We believe this to be a strength in our case since both shape and size of ventral branches probably contribute to the lock-and-key mechanism; for example, spiny but short D. yakuba ventral branches may not harm D. santomea females (Kamimura, 2012;Kamimura & Mitsumoto, 2012b).

| Effect of temperature on size and shape
In most insects and other ectotherms, adult body size typically increases with lower temperatures (Angilletta et al., 2004).
Bergmann's rule, which posits an increasing body size with higher altitude, has been observed within the São Tomé island for the terrestrial caecilian Schistometopum thomense, over a temperature range of 9°C (Measey & Van Dongen, 2006). In contrast to other body parts, the genitalia of insects, and of D. melanogaster in particular, have been reported as not, or little, plastic in response to temperature or other types of environmental variation (Eberhard, 2009;Shingleton et al., 2009;Wheeler et al., 1993). We find here that this is true also for D. yakuba but not for D. santomea: changing the rearing temperature from 25°C to 18°C leads to an increase in spine thrust in D. santomea male genitalia that is similar to what is observed between D. santomea and D. yakuba. In the present study, we only analyzed the effect of temperature on spine thrust, a scalar quantifier which captures one and only one characteristic of the whole shape. We did not examine whether the plasticity-induced change is affecting the entire shape of the ventral branches in the same way as the interspecific change. It would be interesting to include additional landmarks to capture F I G U R E 4 Parameter adjustment for the machine detection algorithm. Training the algorithm relies on two layers of transformation that are each dependent on one parameter: contour coordinate smoothing (horizontal axis) and curvature profile smoothing (vertical axis). Training was performed using a set of 61 individuals, 31 Drosophila santomea 1563 and 30 Drosophila yakuba Oku (see Table 1) for which we manually digitized both landmarks and contours. For each value of the two smoothing parameters, we performed linear regression of spine thrust from manually digitized landmarks against spine thrust derived from automatically digitized landmarks. The colors and values represent the r 2 from that regression. The value used for all detections is contoured in white the entire shape of the genital structure and compare the changes in shape resulting from temperature variation and from interspecific difference (Noble et al., 2019).
The automatically detected landmarks could in principle be used to calculate the centroid size of the anatomical structure, and then, test whether the changes in ventral branches form triggered by temperature reflect heterogeneity in organ size variation within a strain among temperatures, among strains for a given temperature, or even a combination of both.
For each species, we find that strains raised in the same con- Based on our experiments, we cannot fully rule out plasticity in D. yakuba. It is possible that their genital morphology would be altered in external conditions outside of the specific ones that we assayed here. In any case, we find that in our experimental conditions the plasticity of genital form with respect to temperature is higher in D. santomea than D. yakuba.

| Laboratory observations should be complemented by analysis of wild-caught flies
Tests in the laboratory show that D. santomea flies appear to be poorly adapted to high temperatures (Matute et al., 2009). The optimal temperature for larval survival is 21°C for D. santomea and  (Lachaise et al., 2000), which corresponds to temperatures around 18°C or below (https://en.clima te-data.org/, www.world clim.org/bioclim). However, we note that on the southern slopes of the island a few D. santomea flies have also been collected at lower altitudes (650 m) in the dense mist forest near Rio Queijo (Matute & Coyne, 2010;Nagy et al., 2018). This suggests that D. santomea flies can also inhabit warmer regions of the island and that they might be found across the native forest of Sao Tomé, which goes down to sea level on the western slope of the island (Bell & Irian, 2019

| Evolution of the plasticity of ventral branches form
To understand the relevance of this temperature sensitivity of genital form for the past and present evolution of D. santomea and D.
yakuba, more needs to be learnt about their ecology and the plasticity of the ventral branches form of their closely related species, D.
teissieri. Ventral branches are only found in the three species of the D. yakuba complex, D. santomea, D. yakuba, and D. teissieri (Yassin & Orgogozo, 2013). Since ventral branch form plasticity has not been studied in D. teissieri, it is unclear whether this plasticity to F I G U R E 7 Regression tree for spine thrust measures of all Drosophila santomea and Drosophila yakuba isofemale strains at both 18°C and 25°C. Each node gives the spine thrust mean of all samples included in that node and the proportion of the total dataset included in that node. Below each node are two alternatives: to the left the condition is true and to the right the condition is false. Note that the split between D. santomea and D. yakuba happens at the top, thereby suggesting that neither temperature nor years have an effect on spine thrust within D. yakuba

F I G U R E 8 Density distribution of spine thrust values for
Drosophila santomea lines at 18 and 25°C and Drosophila yakuba lines. Density distribution inferred from D. santomea males raised at 18°C (blue) and 25°C (red), D. yakuba males at both temperatures (orange). Distributions are inferred from the total data shown in Figure 2a temperature is an ancestral trait which has been lost in D. yakuba or if it is a novel trait which evolved in D. santomea only. The species D. teissieri is not found in São Tomé but on the mainland and a few islands of the African continent; it can hybridize with D. yakuba (Cooper et al., 2018;Turissini & Matute, 2017). In D. teissieri males, the spines are very long and no layer of cuticle is present between them (Kamimura & Mitsumoto, 2012a;Yassin & Orgogozo, 2013). In any case, even if the extent of ventral branch form plasticity in D.
teissieri was known, it would still be difficult to reconstruct ancestral trait states based on only three species.
The female protective pouches, into which the spiny ventral branches of D. yakuba males fit during copulation, were observed in D. yakuba but not in D. santomea females raised at 21°C and 25°C (Kamimura & Mitsumoto, 2012b;Yassin & Orgogozo, 2013 Such a scenario is opposite to the most common view that posits that morphological diversification tends to proceed through losses of plasticity, rather than gains of plasticity ("flexible stem hypothesis"; Schneider & Meyer, 2017;West-Eberhard, 2003; "plasticityfirst" model; Levis & Pfennig, 2016). It is possible that the decrease in spine thrust that occurred during evolution in the lineage leading to D. santomea was accompanied by a gain of ventral branches form plasticity toward temperature. It is unclear whether the plasticity of ventral branches form to temperature is adaptive. More knowledge about the ecology of D. santomea and its sister species will be required to elaborate a convincing scenario to interpret the role of the ventral branch form plasticity that we discovered.

| CON CLUS ION
Our data show that genitalia can be plastic to temperature and that this plasticity can evolve coincidentally with speciation. Whereas the sensitivity of insect genitalia shape to temperature or nutrition has been used previously as a proof against the lock-and-key hypothesis

ACK N OWLED G M ENTS
We are grateful to São Tomé authorities for allowing us to collect flies. We thank David Stern and Daniel Matute for fly strains. We thank the Courtier laboratory for helpful discussions. We also thank Philippe Rinaudo for feedback on the statistical analyses. The research leading to this paper has received funding from the European

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.