Pigmentation in Drosophila melanogaster reaches its maximum in Ethiopia and correlates most strongly with ultra-violet radiation in sub-Saharan Africa

Fly populations

Table 1 shows the list of the 30 geographical populations that we used in our study. These populations were classified into two major clans: sub-Saharan and cosmopolitan (26 and 4 populations, respectively). The sub-Saharan clan was further divided into four subclans according to their population genetic structure revealed from recent genomic data [30]: west, east, ethiopian and south African subclans (11, 5, 4 and 6 populations, respectively; Figure 1).

Population	Symbol	Latitude	Longitude	Altitude	Date	Data set	N
Sub-Saharan Africa
West
Dondé, Guinea	GU	10.70	−12.25	801	06/2004	I	21
Cotonou, Benin	BN	6.35	2.43	52	05/2004	I	8
Maiduguri, Nigeria	NG	11.85	13.16	295	09/2004	I	9
Kareygorou, Niger	NR	13.55	2.03	195	12/2004	I	20
Yokadouma, Cameroon	CY	3.52	15.05	561	04/2004	I	8
Mbalang-Djalingo, Cameroon	CD	7.32	13.73	1213	03/2004	I, D	20, 6
Maroua, Cameroon	CM	10.60	14.32	385	03/2004	I	22
Nkouondja, Cameroon	CN	5.50	10.68	1121	04/2004	I, D	24, 6
Mbengwi, Cameroon	CW	6.02	10.00	1274	04/2004	I, D	6, 4
Oku, Cameroon	CO	6.25	10.43	2169	04/2004	I, D, M	9, 5, 10
Franceville, Gabon	GA	−1.65	13.60	332	03/2002	I	19
East
Namulonge, Uganda	UG	0.53	32.60	1134	04/2005	I, D, M	20,21,17
Cyangugu, Rwanda	RC	−2.29	28.55	1602	12/2008	D	14
Gikongoro, Rwanda	RG	−2.49	28.92	1927	12/2008	D, M	25, 25
Marigat, Kenya	KR	0.47	35.98	1062	01/2009	D	25
Malindi, Kenya	KM	−1.43	40.03	78	01/2009	D	21
Ethiopia
Gambella, Ethiopia	EA	8.25	34.59	525	12/2011	M	6
Ziway, Ethiopia	EZ	7.93	38.72	1642	12/2008	D	25
Dodola, Ethiopia	ED	6.98	39.18	2492	12/2008	D	7
Fiche, Ethiopia	EF	9.81	38.63	3070	12/2011	M	44
South
Mwanza, Malawi	MW	−15.62	34.52	618	07/2001	I, D	13, 11
Siavonga, Zambia	ZI	−16.54	28.72	530	07/2010	M	51
Lake Kariba, Zimbabwe	ZK	−16.52	28.80	619	05/1994	I	16
Sengwa, Zimbabwe	ZS	−18.16	28.22	865	09/1990	I	13
Phalaborwa, South Africa	SP	−23.94	31.14	375	07/2010	M	6
Dullstroom, South Africa	SD	−25.42	30.10	2000	12/2011	M	49
Cosmopolitan
Palearctic
Cairo, Egypt	EG	30.10	31.32	25	01/2011	M	9
Lyon, France	FR	45.77	4.86	175	07/2010	M	71
Neotropical
Puerto Maldonado, Peru	PM	−12.60	−69.19	200	04/2013	M	14
Cusco, Peru	PC	−13.52	−71.97	3400	04/2013	M	10

N = number of lines.

Populations were also divided into three sets of data according to the location and date of pigmentation scoring: I (Cornell University, Ithaca 2005), D (University of California, Davis 2009) and M (University of Wisconsin, Madison 2013). The I data set includes populations whose abdominal pigmentation scores have already been published by Pool and Aquadro [13]. Among the 19 sub-Saharan populations studied by these authors, we excluded three populations from Eritrea, Kenya and South Africa due to elevated genomic evidence of recent admixture from cosmopolitan strains [30]. The I and D data sets consisted of 228 and 170 isofemale lines, whereas M consisted of 312 lines that were inbred for eight generations (Table 1).

Scoring pigmentation

For each line, five males and ten females were maintained at 20°C and ~75% humidity on standard Drosophila medium (containing molasses, corn meal, yeast, agar, and antimicrobial agents). One, 3–5 day old female from the progeny per line was photographed using an Amscope SM-4TZZ-144A dissection microscope under CO₂ anesthesia. For the I and D data sets flies were photographed on lateral view, but for M both lateral and dorsal views were taken. Photos were then analyzed using the ImageJ software package [42].

We scored different pigmentation traits on two cuticular structures: the fourth abdominal segment and the (thoracic) mesonotum (Figure 2). For abdominal pigmentation two traits were measured: pigmentation intensity near the anterior margin (A4 background), measured in greyscale % as in Pool and Aquadro [13], and width of the posterior black stripe (A4 stripe) as in David et al. [12]. For the latter trait, we divided the width in pixels of the black stripe by the whole width in pixels of A4, rather than visually estimating the ratio. Both traits were measured approximately half-way between the dorsal midline and the dorsal-ventral (tergite-sternite) boundary. A4 stripe was measured for flies from the M data set only. For the thorax, three pigmentation intensity traits were scored in greyscale %: on the latero-anterior margin of the scutum near the humerus, on the posterior margin of the scutum near the scuto-scutellar suture, and on the mesopleuron on the katepisternal sclerite above the level of the two large katerpisternal bristles (hereafter MPL). In addition, we scored the thoracic trident according to the four phenotypic classes proposed by David et al. [14]: 0 = no trident, 1 = faint trident, 2 = clearly marked trident, and 4 = dark trident. Because humeral, presutural and trident pigmentation traits need to be scored from a dorsal perspective they were only scored for flies from the M data set, whereas MPL was scored for all data sets.

Exposure time, zoom width and illumination level using an Amscope adaptor for LED lamp at maximum lighting were kept constant within each data set but they differed due to different laboratory equipments between the three sets. Flies were microphotographed on a fixed region of the CO₂ pad, to minimize background effects. Since some strains were measured in different sets (Table 1), their A4 background means were used to scale greyscale measurements of D and M relative to I using a linear model. To scale D relative to I, three populations were used (MW = 53.6 and 62.6, UG = 68.2 and 74.0, and RG = 61.4 and 72.1, with scores indicating population means for measurements made in D and I, respectively). This led to a relationship of x _I = 21.23 + 0.79 x _D (R ² = 0.89), with x _D being the mean of measurements scored in Davis and x _I the scaled measurement relative to Ithaca. To scale M relative to I, three data points were also used, corresponding to CO, UG and ZK/ZI, with the latter being two populations from Zimbabwe and Zambia, respectively, that are less than 10 km apart and share similar habitats. The pigmentation scores (CO = 62.0 and 72.1, UG = 59.4 and 74.0, and ZI/ZK = 41.8 and 61.8, with measurements indicating means for M and I, respectively) led to a relationship of x _I = 37.92 + 0.58 x _M (R ² = 0.93), with x _M being the mean of measurements scored in Madison and x _I the scaled measurement relative to Ithaca. These rescalings are unlikely to perfectly account for the differences between data sets, which may add noise to some analyses. However, we do not expect any strong bias in terms of environmental correlations, and we find no significant effect of data set on rescaled A4 and MPL pigmentation scores for rescaled A4 and MPL pigmentation scores, respectively: Kruskal-Wallis chi-squared = 4.46 (P = 0.1073) and 4.19 (P = 0.1233).

Statistical analyses

Three main statistical analyses were conducted in this paper. First, we evaluated the levels of phenotypic correlations between the five pigmentation traits scored within two sub-Saharan populations representing phenotypic extremes of the M data set. To correct for multiple correlation tests, we estimated the false discovery rate (FDR) for a given P value, i.e. a q value, and significance levels were considered at q < 0.05, i.e. a 5% FDR cutoff.

where y is the pigmentation score mean of a population, a is the intercept, b ₁ and b ₂ are the slopes (or the clines) of the absolute latitude and altitude values, and ε is the residual. Analyses were conducted for sub-Saharan populations only. For comparative purposes, we also independently estimated the zeroth-ordered correlation coefficient (r₀) between each pigmentation trait, latitude and altitude, as well as the semi-partial correlation coefficient (r_s) of each pigmentation trait with the latitude or the altitude after controlling for the effect of the other geographical variable on pigmentation.

Third, we analyzed the correlation between A4 background and MPL and different ecological and historical factors that may affect pigmentation evolution. In order to determine the environmental factors that may explain each of the clines, we obtained GIS meteorological and geological data for each population. Meteorological data averaged over 22 years (from 1983 to 2005) were extracted as annual averages from the NASA Surface meteorology and Solar Energy: Global Data Sets website (www.eosweb.larc.nasa.gov). These included five radiation sets (namely, insolation incident, diffuse radiation, direct normal radiation and latitude tilt radiation, each measured in kWh/m²/day) and seven climatic sets (namely, average, minimum and maximum air temperatures and earth temperature measured in °C, relative humidity in %, atmospheric pressure in kPa, and wind speed at 50 m above ground in m/s). We also included annual average values of UV index from the Tropospheric Emission Monitoring Internet Service (www.temis.nl) in units of 25 mW/m². Geological data of soil sand, silt and clay contents at <2 μm in % were obtained from the Africa Soil Information Service (www.ciesin.columbia.edu/afsis). We estimated the coefficients of correlation between pigmentation traits and each variable, and then chose for each ecological category that has previously been hypothesized to affect pigmentation (i.e. UV protection, thermal budget, desiccation and crypsis) a single variable showing the strongest correlation with pigmentation. In addition we included two principal component scores, namely PC1 and PC2, from Pool et al. [30] reflecting subclan differentiation to control for population structure and historical effects. To evaluate the relative importance of each of these ecological or historical factors, we used an information-theoretic approach that ranks a number of linear model-based hypotheses according to evidence ratios and posterior probabilities [43]. We conducted this analysis first using each factor as a separate model, and then using each factor and combination of two factors as models. For each model, the sample-size corrected Akaike information criterion (AICc) was estimated, and models were ranked on the basis of increasing AICc and evaluated in respect to the descending probabilities (w).

All statistical analyses were performed using the R software package (www.r-project.org). Semi-partial and second-ordered partial correlation analyses were conducted using the ppcor 1.0 package [44] in R. False Discovery Rate (FDR) control [45] was applied for multiple tests of phenotypic correlations between traits using the LBE 1.22 software package [46] in R. For model set evaluation, we modified the R code provided by Correa and Hendry [47] to fit our data.

Correlation between pigmentation traits

A new phenotype of thoracic pigmentation in Ethiopia

A uniquely melanic phenotypic class which goes beyond the scale defined by David et al. [14] was found and scored for some lines from the high-altitude Ethiopian population EF. These phenotypes have a diffuse wide trident on the thorax making the mesonotum (and the head) appearing completely black (Figure 2). They were scored as a fifth class (value = 4), a phenotype that was absent in any previously examined sub-Saharan or cosmopolitan strain of D. melanogaster or its close relative D. simulans (J. R. David, pers. comm.). Consequently, the EF mean of thoracic trident was nearly twice that of other high latitude or altitude populations (e.g., France, Peru, Cameroon, Uganda and Rwanda) (Figure 3). The phenotype was also found in another high altitude Ethiopian population (ED), but it was absent from the low-altitude Ethiopian population EA.

Geographical differentiation of pigmentation traits

A4 background and MPL pigmentation variation between populations in all data sets after calibration with the scale of the I data set showed similar geographical patterns between genetic clans (ANOVA's F = 3.09 (P = 0.0134) vs. 3.22 (P = 0.0292), respectively; Figure 4). The two traits did not differ, however, in their variance, with the coefficient of variation being 6.74% for A4 background and 5.63% for MPL (F-test = 1.43; P = 0.83, [48]). Tukey's pairwise test showed a single significant pigmentation difference between Ethiopian and south African populations (Table 3; P = 0.0127 and 0.0152 for A4 background and MPL, respectively). Remarkably, the uniquely melanic thoracic phenotype of high-altitude Ethiopian populations found at altitudes ranging between 1,600 and 3,100 m was not found in any other high-altitude sub-Saharan or even cosmopolitan population, including the population from Cusco, Peru (3,400 m). Melanism from Ethiopian highland strains thus appears to greatly exceed the pigmentation of any previously described D. melanogaster populations.

Opposite latitudinal clines between sub-Saharan and cosmopolitan populations

A4 background and MPL pigmentation in sub-Saharan Africa were significantly correlated with both latitude and altitude (Table 4; Figure 5). This situation is similar to what has previously been found on the Indian subcontinent [15],[19]. Although in both cases pigmentations' slopes were positive for altitude, i.e. populations inhabiting high altitudes are darker, a negative slope was found for latitude in sub-Saharan Africa, whereas the latitudinal slope was positive in India. This difference is most likely due to the topographic nature of the Indian subcontinent where both latitudes and altitudes covary, i.e. elevations increase northward toward the Himalayan plateau. In Africa, latitudes and altitudes do not covary (r = −0.12, P = 0.56) and the pigmentation clines with each factor persist after controlling for the other factor (r_latitude = −0.64 (P = 6.7 × 10⁻⁵) and −0.52 (P = 0.004), r_altitude = 0.75 (P = 7.2 × 10⁻⁸) and 0.59 (P = 4.8 × 10⁻⁴), for A4 background and MPL, respectively).

Variable	Intercept	Latitude	Altitude
A4 background	65.866 ± 1.209***	−0.321 ± 0.081***	0.004 ± 0.001***
MPL	71.859 ± 1.336***	−0.257 ± 0.089**	0.003 ± 0.001**

*P < 0.05; **P < 0.01; ***P < 0.001.

Pigmentation correlates most strongly with UV radiation in sub-Saharan Africa

Among the 18 environmental and historical factors considered in this study, pigmentation significantly correlated with 18 and 11 variables for A4 background and MPL, respectively (Table 5; Figure 6). For both traits UV index had the strongest correlations (r = 0.77 and 0.78 for A4 and MPL, respectively; Figure 6A). In order to evaluate the different ecological and historical hypotheses, we conducted a model rank analysis using five environmental and two population genetic structure factors. We found UV radiation to be the most explanatory factor for variation in sub-Saharan Africa, when each environmental factor is considered separately (g₀₁ in Tables 6 and 7). Flies are darker at higher levels of UV index (slope = 3.95 ± 0.74 (P = 1.7 × 10⁻⁵) and 3.58 ± 0.62 (P = 6.6 × 10⁻⁶), for A4 background and MPL, respectively). Soil silt content (g₀₅) was always the second most explanatory factor (slope = 0.35 ± 0.07 (P = 4.3 × 10⁻⁵) and 0.28 ± 0.07 (P = 3.1 × 10⁻⁴), for A4 background and MPL, respectively), but for both traits it was less likely than UV index (evidence ratio: E_g01,g05 = 2.6 and 49.0, for A4 background and MPL, respectively). Similar results were obtained when combinations of traits were analyzed in a multivariate regression context (Tables 8 and 9). For MPL, the first four models, which weighted for ~75% of the variance among models, were models including UV index alone or in combinations (Table 8). For A4 background, the most explanatory model which weighted for 61%, did not include UV index. Instead this model was explained by correlations with temperature and the first principal component of genetic differentiation (slopes = −1.08 ± 0.17 (P = 1.97 × 10⁻⁶) and 27.08 ± 5.51 (P = 5.7 × 10⁻⁵), respectively), although the second and third most likely combinations did include UV. Thus, the unique melanic phenotype in Ethiopia may have evolved to protect against UV radiation, which reaches its maximum levels (>300 mW/m²) in the Ethiopian highlands. However, this hypothesis requires further study, as does the potential role of other ecological factors.

	A4 background	MPL
Radiation
UV index	0.737***	0.760***
Insolation incident	−0.157*	−0.005
Diffuse radiation	0.411**	0.300
Direct normal radiation	−0.319*	−0.186
Latitude tilt radiation	−0.279*	−0.136
Clear sky radiation	0.324*	0.389*
Temperature
Average temperature	−0.596***	−0.495**
Average minimum temperature	−0.485**	−0.413*
Average maximum temperature	−0.640***	−0.532**
Earth temperature	−0.665***	−0.501**
Humidity
Relative humidity	0.414**	0.256
Atmospheric pressure	−0.658***	−0.607**
Wind
Wind speed	−0.305*	−0.227
Soil
Sand	−0.612***	−0.520**
Silt	0.714***	0.652**
Clay	0.178*	0.113
Genetic structure
PC1	0.465**	0.491**
PC2	−0.323*	−0.376*

Significance levels after FDR correction.

*q < 0.05; **q < 0.01; ***q < 0.001.

Pigmentation as a modular trait

Insect pigmentation can be regarded as a modular trait: its development depends on a conserved network of structural genes of the melanin synthesis pathway whose expression levels may differ between different body organs [38]. We investigated this hypothesis within two segments: the 2^nd thoracic segment (mesonotum) and the 4^th abdominal segments (A4), which are developmentally four segments apart. Previous studies in Drosophila showed that the degree of genetic correlation of pigmentation decreases with increasing distance between body segments [49],[50]. In agreement with these studies, we did not find a significant correlation between the two segments in the lowland ZI population which is presumably under no directional selection for pigmentation. On the other hand, pigmentation traits correlated between the two segments in the high-altitude Ethiopian population, potentially due to selectively favored alleles that alter multiple pigmentation traits. A recent experimental selection study has shown the presence of correlated response to selection between the thorax and the 2^nd abdominal segment [41]. Geographical variation of A4 pigmentation in sub-Saharan Africa and of thoracic trident in Australia were found to be related to differences in ebony expression levels [7],[13],[16], but detailed molecular dissection of the regulation of ebony expression revealed that different enhancers of the gene affect different body parts [7],[10],[51].

Adaptive significance of Drosophilamelanism

Latitudinal or altitudinal clines of abdominal and thoracic pigmentation has previously been found in D. melanogaster[13]–[20] and also in other drosophilid species [21]–[25]. In all cases darker phenotypes were encountered at high latitudes or altitudes, but the exact cause of the clinal variation remained elusive. The most invoked hypothesis was the thermal budget or thermal melanism stating that darker flies absorbs better solar radiation in colder environments [32]. This might explain some longitudinal clines too [52]. However, the high surface to volume ratio of these small insects may preclude them from maintaining a higher body temperature than their surroundings [53]. There is little experimental evidence that dark flies are warmer [54], and no experiment to our knowledge has shown a higher fitness of dark morphs in colder environments. The second most cited hypothesis is desiccation resistance. In India, variation in desiccation resistance parallels latitudinal and altitudinal clines of pigmentation [55], and there is experimental evidence that dark morphs are more resistant to desiccation than light morphs due to cuticular thickness [19]. However, the altitudinal cline of pigmentation in sub-Saharan D. yakuba was not associated with desiccation resistance [25].

We conducted intensive correlation studies between clinal variation in pigmentation and different environmental factors and found that UV radiation was the strongest predictor of pigmentation traits especially for thoracic pigmentation (MPL). For abdominal pigmentation (A4 background), UV was also the strongest explanatory factor when factors were analyzed separately, but combinations of other factors such as temperature and population history might also be relevant. We did not find any relevance to atmospheric pressure or relative humidity, the two factors which may be the most explanatory for a desiccation-resistance hypothesis. In endotherms, pigmentation tends to increase with decreasing latitudes, the so-called Gloger's rule [5] (cf. [56] for an example in humans). In Drosophila, the opposite was always encountered, i.e. a positive relationship between pigmentation and latitude, in Europe, India and Australia [14]–[16],[18]. However, David et al. [14] noted that latitudinal clines for thoracic pigmentation were only found in the temperate region encompassing North Africa and Europe (above 30°), whereas pigmentation in India also correlated with altitude which covaries with latitude [15]. Pool and Aquadro [13] found a negative latitudinal cline in sub-Saharan Africa, but the cline was not significant when they corrected for altitudinal effect. In our study, the negative latitudinal cline persisted even after correcting for altitudes. This may be due to our inclusion of high-altitude populations from both higher and lower latitudes. Our study thus provides the first evidence of a negative latitudinal cline of abdominal and thoracic pigmentation in sub-Saharan Africa in agreement with Gloger's rule but for ectotherms.

Considering altitudinal clines of pigmentation, it has long been suggested that the “ecological importance of melanism at high altitudes is closely bound up with the fact that dense pigments serve as effective protection against the injurious effects of the intense ultra-violet” (p. 53, [57]). Scott [58] was “struck by the unicolorous black hue of the body” of some snails and one chloropid fly in high altitudes in Ethiopia, and we report here a uniquely dark morph of D. melanogaster from Ethiopia as well. Pool and Aquadro [13] demonstrated the presence of altitudinal clines of abdominal pigmentation in sub-Saharan Africa. Curiously, there is no obvious altitude-pigmentation cline in D. simulans (JEP, personal observation), even though this species is thought to spend more time outdoors than its close relative [59], and exists at high altitudes in Africa (e.g., Dodola, Ethiopia, where pigmentation is essentially a diagnostic trait between D. simulans and D. melanogaster).

Laboratory exposure of D. melanogaster to elevated UV radiation has typically shown wild type flies to be more resistant than light or dark pigmentation mutant strains [34]–[36], but flies experimentally fed melanin were found to be more resistant to irradiation [60]. Recently, Matute and Harris [25] reported the association between an altitudinal cline of pigmentation in Drosophila yakuba and UV resistance on two islands near western Africa. Surprisingly, the authors found that lightly pigmented strains had the longest survival after exposure to high levels of UV radiation. In contrast, the correlations we observe are consistent with a protective effect of dark pigmentation against UV radiation. These species could conceivably differ in the relationship between pigmentation and UV protection. Another possibility is that resistance to elevated UV doses in the laboratory may differ in important ways from UV resistance under natural conditions (e.g. due to the release of heat by melanin exposed to high levels of UV in the lab). Future UV resistance experiments with light and dark populations of D. melanogaster should be conducted to explain clinal variation in sub-Saharan Africa.

Our analyses were conducted on flies raised under similar laboratory conditions, and hence investigate genetic differences between populations. However, the phenotypes that occur in nature may be modulated by phenotypic plasticity [61], such as the documented influence of developmental temperature on D. melanogaster pigmentation [22],[49],[62]–[68], which could also differ between populations. Hence, it will also be worthwhile to study the influence of temperature and other environmental factors on the pigmentation of D. melanogaster from Ethiopia and elsewhere.