GitHub respository for manuscript entitled ‘Continent-wide differentiation of fitness traits and patterns of climate adaptation among European populations of Drosophila melanogaster’

Author contributions can be found here

Contents of the repository, raw data, tables, figures and scripts can be found here

Correspondence: esra@bmb.sdu.dk; thomas.flatt@unifr.ch

Short Running Title: Phenomic Analysis of European D. melanogaster

Lay Summary

For over 100 years, the vinegar fly Drosophila melanogaster has been a workhorse for studying the genetic and phenotypic basis of evolutionary change. Ancestrally of southern-central African origin, this insect expanded its range, migrated out of Africa and colonized new habitats and climate zones around the globe. Recent population genomic analyses suggest that Europe was colonized from the Middle East ~1,800 years ago, yet we still know surprisingly little about how this human commensal adapted to different locales on the European continent and elsewhere. Here, we assayed fly strains isolated from 9 European populations for 16 traits in an international consortium effort, involving >100 researchers from 26 research groups in 17 countries and using semi-standardized experimental protocols. Despite differences in experimental (environmental) conditions among labs, our trait measurements agreed well across research teams, suggesting that they are robust and generalizable. European populations of the vinegar fly are markedly differentiated for several phenotypic traits that are related to Darwinian fitness and which might be subject to natural selection. Notably, flies from locations with higher humidity and rainfall and lower temperature survive better to the adult stage, lay more eggs, are better at surviving famine and, as adults, live longer than flies from drier, warmer locations; yet, as compared to flies from warmer, drier locales, they are considerably less well adapted to surviving heat stress, suggesting a trade-off between these components of fitness. Our findings begin to illuminate how this subtropical/tropical insect has adapted to different climates and habitats on the European continent and showcase the value of collaborative multi-lab studies that perform experiments in parallel.

1 - Materials and Methods

1.1 General overview of the phenotyping effort

We assayed 16 phenotypic traits, most of them representing fitness components (see Table 1 in main text for details): viability, developmental time, dry weight, thorax length, wing area, fertility, lifespan, cold-shock mortality, chill-coma recovery time, heat-shock mortality, diapause, locomotor activity, circadian eclosion timing, pigmentation, starvation resistance, and parasitoid resistance. Assay protocols for each trait are given below; detailed information about phenotyping batches, experimental blocks, replication, sample sizes, etc. can be found here https://github.com/esradm/DrosEU_PhenotypingWG/tree/main/SummaryTables. Phenotyping assays were carried out at 25ºC, 12 hours light:12 hours dark, and a minimum relative air humidity of 60%. Strains were maintained under density-controlled conditions for at least two generations prior to the assays. The fly food recipes used by the different labs are shown in section 1.3. Unless stated otherwise, mated females and males were phenotyped for a given trait. The following principal investigators and their research groups participated in the phenotyping effort (in alphabetical order): Jessica Abbott (JA); Alan Bergland (AB); Jean-Christophe Billeter (JCB); Hervé Colinet (HC); Claudia Fricke (CF); Thomas Flatt (TF); Patricia Gibert (PG); Josefa González (JG); Sonja Grath (SG); Katja Hoedjes (KH); Jan Hrcek (JH); Iryna Kozeretska (IK); Julian Mensch (JM); Banu Onder (BO); John Parsch (JP); Elena Pasyukova (EP); Nico Posnien (NP); Michael G. Ritchie (MR); Christian Schlötterer (CS); Paul Schmidt (PS); Marina Stamenkovic-Radak (MSR); Eran Tauber (ET); Jorge Vieira (JV); Christian Wegener (CW); Bas J. Zwaan (BZ).

Nine sampling locations were chosen based on genomic data and these locations covered a wide range of latitude (~20°) and longitude (~40°) across the continent (see Table S1 below). From each location, 15 to 20 isofemale lines were established in corresponding labs at the sampling location and the isofemale lines were centrally maintained by Élio Sucena at Instituto Gulbenkian de Ciência (IGC), Lisbon, Portugal. A total of 173 isofemale lines were used in this study.

Table S1: Sampling locations, coordinates and collectors, sorted by increasing latitude. Altitude is expressed in meters.

Country	Location	LocationAbbr	Latitude	Longitude	Altitude	Collector	CollectionDate	DEST2.0_Name	DEST2.0_locality_Code
Portugal	Recarei	RE	41,15	-8,41	175	Jorge Vieira	05/10/2018	PT_Por_Rec_1_2018_10_05	PT_Por_Rec
Spain	Gimenells (Lleida)	GI	41,618	0,62	173	Josefa Gonzalez, Marta Pascual	25/08/2018	ES_Ler_Gim_1_2018_08_25	ES_Ler_Gim
Denmark	Karensminde	KA	55,945	10,213	15	Mads Schou	05/10/2018	DK_Mid_Kar_1_2018_10_05	DK_Mid_Kar
Germany	Munich	MU	48,18	11,61	520	John Parsch, Amanda Glaser-Schmitt	21-30/06/2018	DE_Bay_Mun_1_2018_06_25	DE_Bay_Mun
Austria	Mauternbach	MA	48,375	15,56	572	Andrea Betancourt	08/09/2018	AT_Nie_Mau_1_2018_09_08	AT_Nie_Mau
Finland	Akaa	AK	61,1	23,52	88	Maaria Kankare	20/07/2018	FI_Pir_Aka_1_2018_07_20	FI_Pir_Aka
Ukraine	Uman	UM	48,753	30,206	214	Iryna Kozeretska	18/08/2018	UA_Che_Uma_1_2018_08_18	UA_Che_Uma
Turkey	Yesiloz	YE	40,231	32,26	680	Banu Onder	27/09/2018	TR_Ank_Yes_1_2018_09_27	TR_Ank_Yes
Russia	Valday	VA	57,979	33,244	217	Elena Pasyukova	20-30/08/2018	RU_Nov_Val_1_2018_08_25	RU_Nov_Val

Table S2. Contributing labs and traits, sorted by trait. Author contributions can be found here

Country	City	Supervisor.PI	Trait	Contributors
Argentina	Buenos Aiers	Mensch	Chill-coma recovery time	Florencia Putero, Lucas Kreiman, Julian Mensh
Portugal	Porto	Vieira	Chill-coma recovery time	Jorge Vieira, Cristina P. Vieira, Pedro Duque, Tânia Dias
Germany	Würzburg	Wegener	Circadian eclosion timing	Susanne Klühspies, Christian Wegener
Spain	Barcelona	Gonzalez	Cold-shock mortality	Llewellyn Green, Josefa Gonzalez, Miriam Merenciano
Ukraine	Kyiv	Kozeretska	Cold-shock mortality	Svitlana Serga, Alexandra Protsenko, Oleksandr Maistrenko, Iryna Kozeretska
Portugal	Porto	Vieira	Cold-shock mortality	Jorge Vieira, Cristina P. Vieira, Pedro Duque, Tânia Dias
France	Lyon	Gibert	Development time	Cristina Vieira, Laurence Mouton, Natacha Kremer, Sonia Martinez, Patricia Gibert
Germany	Munich	Grath	Development time	Ingo Müller, Sonja Grath
The Netherlands	Lausanne	Hoedjes	Development time	Hristina Kostic, Katja Hoedjes
USA	Philadelphia	Schmidt	Development time	Ozan Kiratli, Yonatan Babore, Liam Forsythe, Paul Schmidt
Serbia	Belgrade	Stamenkovic-Radak	Development time	Marija Savic Veselinovic, Marija Tanaskovic, Aleksandra Patenkovic, Mihailo Jelic, Katarina Eric, Pavle Eric, Slobodan Davidovic, Marina Stamenkovic-Radak
The Netherlands	Wageningen	Zwaan	Development time	Joost van den Heuvel, Bas Zwaan
USA	Charlottesville	Bergland	Diapause	Liam Miller, Alan Bergland, Priscilla Erickson
Switzerland	Fribourg	Flatt	Diapause	Esra Durmaz, Envel Kerdaffrec, Thibault Schowing, Virginie Thieu, Marisa Rodrigues, Thomas Flatt
Austria	Vienna	Schlötterer	Diapause	Manolis Lyrakis, Christian Schlötterer
France	Rennes	Colinet	Dry weight	Sapho-Lou Marti , Hervé Colinet
The Netherlands	Lausanne	Hoedjes	Dry weight	Hristina Kostic, Katja Hoedjes
Turkey	Ankara	Onder	Dry weight	Seda Coskun, Senel Selin Senkal, Dogus Can, Banu Sebnem Onder
The Netherlands	Groningen	Billeter	Fertility	Xiaocui Wang, Tiphaine Bailly, Mario Mira, Jean-Christophe Billeter
Germany	Muenster	Fricke	Fertility	Claudia Fricke
Portugal	Lisbon	Sucena	Fly husbandry	Tânia Paulo, Elio Sucena
Germany	Munich	Parsch	Heat-shock mortality	Eliza Argyridou, Amanda Glaser-Schmitt, John Parsch
Portugal	Porto	Vieira	Heat-shock mortality	Jorge Vieira, Cristina P. Vieira, Pedro Duque, Tânia Dias
Switzerland	Fribourg	Flatt	Lifespan	Esra Durmaz, Envel Kerdaffrec, Thibault Schowing, Virginie Thieu, Marisa Rodrigues, Thomas Flatt
Germany	Munich	Parsch	Lifespan	Amanda Glaser-Schmitt, Eliza Argyridou, John Parsch
Russia	Moscow	Pasyukova	Lifespan	Natalia Roshina, Alexander Symonenko, Mikhail Trostnikov, Evgenia Tsybul’ko, Ekaterina Veselkina, Olga Rybina, Elena Pasyukova
Israel	Haifa	Tauber	Locomotor activity	Bettina Fishman, Eran Tauber
Czech Republic	Ceske Budejovice	Hrcek	Parasitoid resistance	Vincent Montbel, Somayeh Rasouli Dogaheh, Jan Hrcek
Sweden	Lund	Abbott	Pigmentation	Jessica Abbott, Qinyang Li, Shahzad Khan
France	Lyon	Gibert	Pigmentation	Cristina Vieira, Laurence Mouton, Natacha Kremer, Sonia Martinez, Camille Mermet, Patricia Gibert
USA	Philadelphia	Schmidt	Pigmentation	Amy Goldfischer, Paul Schmidt
Spain	Barcelona	Gonzalez	Starvation resistance	Llewellyn Green, Josefa Gonzalez, Miriam Merenciano
Turkey	Ankara	Onder	Starvation resistance	Seda Coskun, Ekin Demir, Senel Selin Senkal, Cansu Aksoy, Banu Onder
Russia	Moscow	Pasyukova	Starvation resistance	Alexander Symonenko, Natalia Roshina, Mikhail Trostnokov, Ekaterina Veselkina, Evgenia Tsybul’ko, Olga Rybina, Elena Pasyukova
Ukraine	Kyiv	Kozeretska	Thorax length	Svitlana Serga, Alexandra Protsenko, Oleksandr Maistrenko, Iryna Kozeretska
Germany	Göttingen	Posnien	Thorax length	Micael Reis, Lennart Hüper, Nico Posnien
UK	St Andrews	Ritchie	Thorax length	Megan Mcgunnigle, Nicola Cook, Teresa Abaurrea, Michael Ritchie
USA	Philadelphia	Schmidt	Thorax length	Amy Goldfischer, Paul Schmidt
USA	Philadelphia	Schmidt	Time to pupation	Paul Schmidt
France	Lyon	Gibert	Viability	Cristina Vieira, Laurence Mouton, Natacha Kremer, Sonia Martinez, Patricia Gibert
Germany	Munich	Grath	Viability	Ingo Müller, Sonja Grath
The Netherlands	Lausanne	Hoedjes	Viability	Hristina Kostic, Katja Hoedjes
USA	Philadelphia	Schmidt	Viability	Ozan Kiratli, Yonatan Babore, Liam Forsythe, Paul Schmidt
Serbia	Belgrade	Stamenkovic-Radak	Viability	Marija Savic Veselinovic, Marija Tanaskovic, Aleksandra Patenkovic, Mihailo Jelic, Katarina Eric, Pavle Eric, Slobodan Davidovic, Marina Stamenkovic-Radak
The Netherlands	Wageningen	Zwaan	Viability	Joost van den Heuvel, Bas Zwaan
Turkey	Ankara	Onder	Wing area	Cansu Aksoy, Ekin Demir, Ezgi Cobanoglu, Banu Sebnem Onder
Germany	Göttingen	Posnien	Wing area	Micael Reis, Lennart Hüper, Nico Posnien
UK	St Andrews	Ritchie	Wing area	Megan Mcgunnigle, Nicola Cook, Teresa Abaurrea, Marija Tanaskovic, Michael Ritchie
Serbia	Belgrade	Stamenkovic-Radak	Wing area	Marija Savic Veselinovic, Marija Tanaskovic, Aleksandra Patenkovic, Filip Filopovski, Mihailo Jelic, Katarina Eric, Pavle Eric, Slobodan Davidovic, Marina Stamenkovic-Radak

1.2 - Drosophila population samples, isofemale lines, and stock maintenance

In summer and fall 2018 we collected inseminated D. melanogaster females at 9 European locations, which had already previously been characterized at the population genomic level by the DrosEU consortium (see Figure 1 in the main text; see Kapun et al., 2020, 2021; also see Machado et al. 2021): PT, Portugal (Recarei = RE); ES, Spain (Gimenells = GI [Lleida]); TR, Turkey (Yesiloz = YE); DE, Germany (Munich = MU); AT, Austria (Mauternbach = MA); UA, Ukraine (Uman = UM); DK, Denmark (Karensminde = KA); FI, Finland (Akaa = AK); RU, Russia (Valday = VA). From these females we established 173 isofemale lines, with each population being represented by at least 15 lines – we call this collection of lines the DrosEU population panel (DPP). Lines were centralized and maintained in the lab of Élio Sucena (ES) where species status (D. melanogaster vs. D. simulans) was verified using the protocol of Faria & Sucena (2017) and where lines were tested for Wolbachia infection (see below for methodological details). Subsequently, lines were sent to participating laboratories for phenotyping (Figure 1 in the main text and Table S2). To comply with the Nagoya protocol, material transfer agreements (MTAs) were prepared and exchanged between researchers to transport fly samples across borders.

1.3 - Diets

1.4 - Chromosomal inversion frequencies

For each isofemale line, the presence or absence of 5 cosmopolitan inversions (In(2L)t, In(2R)NS, In(3L)P,In(3R)P and In(3R)mo) was diagnosed by PCR using previously published inverted-specific primers and conditions (Corbett-Detig et al, 2012) on DNA extracted from pools of 10-15 mixed-sex flies per line using the Qiagen DNeasy Blood & Tissue Kit. These data were used to investigate the effect of inversions on trait variation (linear model with inversion presence (1) or absence (0) used as explanatory variable). In addition, average inversion frequencies were derived for each population.

1.5 - Trait-specific assay conditions and measurements

Viability

Gibert, Grath, Hoedjes, Schmidt, Stamenkovic-Radak and Zwaan Labs

Viability (=proportion egg-to-adult survival) was determined in parallel with developmental time (see below). Groups of 3-day to 5-day-old adults (at least 25 pairs) per isofemale line were allowed to lay eggs for at least 2 hours. Yeast was provided to stimulate egg laying. Immediately thereafter 40 eggs were collected and placed into a vial with 5 mL medium, with three replicate vials per isofemale line. The PS lab followed an alternative protocol whereby females were allowed to lay eggs directly into vials and the total number of eggs was recorded afterwards. Viability was calculated per vial as the proportion of adult individuals which eclosed from the eggs.

Developmental time

Egg-to-pupa developmental time

Schmidt Lab

Egg-to-pupa developmental time (in hrs) was scored twice per day, when the chamber lights were turned on and two hours before they were turned off, by counting the number of individuals that had pupariated. Locations on the walls of the vials with new pupae were marked with a permanent marker in order to keep track of which pupae had emerged each day.

Egg-to-adult developmental time

Gibert, Grath, Hoedjes, Schmidt, Stamenkovic-Radak and Zwaan Labs

Egg-to-adult developmental time (in hrs) was scored twice per day, similar to egg-to-pupa developmental time (see above); all adults that had emerged from each vial were collected for counting and sexing.

Dry weight

Colinet, Hoedjes and Onder Labs

Measurements were performed in 3, 5 and 5 batches for the HC, KH and BO labs, respectively. Seven-day-old flies (on average 24 females and 24 males per isofemale line) from cultures of controlled density (50 eggs/5 mL of food) were sacrificed by either placing them directly at -20°C or by putting them first into a vial containing a piece of cotton imbibed with an ethyl acetate solution before storing them at -20ºC. Flies were sexed, placed individually into 96-well plates and transferred into a drying oven set at 60-70°C for at least 72 hours. Plates that were not measured immediately after the drying period were stored at room temperature using a protective cover for later weighing; if this was the case, previously dried plates were placed in the oven (60-70°C) for another 24 hours the day before weight measurements, to ensure that flies were well dehydrated. Flies were then transferred individually with an aspirator onto a small piece of aluminum foil for weight measurement on a microbalance with an accuracy of 1µg (Mettler Toledo UMX2 or MT5, Sartorius Cubis Micro Balance).

Thorax length

Kozeretska, Posnien, Ritchie and Schmidt Labs

Flies used for measurements of thorax length were reared in 2 batches by the IK and in 1 batch by all other labs. NP kept the flies at -20°C for 5-14 days prior to mounting. Bodies / thoraces of 5- to 10-day-old flies (per isofemale line and sex, on average ~40 individuals for the IK and MR labs, and ~10 individuals for the NP and PS labs) were placed onto a double-sided sticky tape attached to a microscope slide or taped directly to the slide. Bodies were laid out on their right side, and photographs of thoraces were taken using a digital camera connected to a dissecting microscope (NP lab: QImaging MicroPublisher 5.0 digital camera mounted to a Leica M205 FA stereo microscope; MR lab: Leica DFC295 camera attached to a Leica M60 microscope; IK lab: a digital single-lens reflex (DSLR) Sigeta MCMOS 5100 5.1Mp USB 2.0 camera mounted on a MBS-10 stereo microscope, PS lab: Olympus DP73 digital camera mounted on a Leica MZ9.5 stereo microscope). Within each lab always the same magnification (50x) and resolution were used to increase reproducibility; a scale bar inserted on each photo or a software (PS lab) allowed transforming (square) pixels into (square) µm. Thorax length was defined as the distance from the anterior margin of the thorax to the posterior tip of the scutellum and was measured using the ‘Straight Line’ tool in ImageJ/Fiji (https://fiji.sc/) (Figure S1A).

Wing area

Onder, Posnien, Ritchie and Stamenkovic-Radak Labs

Flies used for wing area measurements were reared in 7, 1, 15 and 5 batches by the BO, NP, MR and MSR labs, respectively. Both the left and right wings of 5- to 10-day-old flies (per isofemale line and sex, on average 30 individuals were measured in the BO, MR and MSR labs, and 10 individuals in the NP lab) were removed and mounted in Entellan (Merck) (BO) or Hoyer’s mounting medium (NP) or placed directly onto double-sided sticky tape (MR, MSR). NP and MSR kept the flies at -20°C for 5 - 14 days prior to wing dissections. Photographs of wings were taken using a digital camera connected to a dissecting microscope (NP: QImaging MicroPublisher 5.0 digital camera mounted to a Leica M205 FA stereo microscope; MSR: Bresser MikroCam 5.0 MP digital camera mounted to Nikon SMZ 745T stereo microscope; MR: Leica DFC295 camera attached to Leica M60 stereo microscope; BO: Leica S9i with integrated 10 MP CMOS-camera). Within each lab the same magnification and resolution was used (NP: 50X; MSR: 50X; MR: 50X; BO: 48X). A photograph of a scale bar (ruler) was taken with the same settings as for the rest of the images during every measurement session, thereby allowing the conversion of (square) pixels to (square) µm.

Figure S1. Thorax and wing measurements. (A) Thorax length was measured from the anterior margin of the thorax to the posterior tip of the scutellum. (B) Measurement area of the wing. (B’) Landmarks used to exclude the posterior part of the wing for area measurements. (C) Position of 15 landmarks used for wing centroid size (WCS) assessment.

Wing area was estimated by quantifying wing centroid size (WCS; Bookstein, 1996). Fifteen landmarks (Figure S1C) were placed on wing images using tpsUtil and tpsDig2 (Rohlf, 2015) to obtain raw x and y landmark coordinates. WCS was calculated as the square root of the sum of squared deviations of landmarks around their centroid as implemented in MorphoJ (Klingenberg, 2011) or tpsRelw (Rohlf, 2015). To confirm that WCS represents wing area, the NP lab also manually measured wing area by outlining wings (Figure S1B) with the “Polygon Selection” tool in Fiji (https://fiji.sc/). The most proximal part of the wing was excluded from measurements, using two distinct landmarks, since it was more susceptible to damage during dissection (Figure S1B’). Linear regression of WCS against manually measured wing area showed a high positive correlation (adjusted R2 = 0.83, F1, 538 = 2605, p < 2.2 x 10-16) (Figure S2).

Figure S2. Linear regression of WCS against manually measured wing area.

Fertility

Billeter and Fricke Labs

For each isofemale line, up to 10 males and 10 virgin females were placed together in single-sex groups and allowed to mature for five days. Individual pairs were then placed together in a vial (5-7 pairs per line), and mating interactions were observed to ensure successful mating (a copulation duration of at least 10 minutes). Up to 5 successfully mated females per isofemale line were obtained. After successful mating, males were discarded and single females were allowed to oviposit for 2 days, then moved to a fresh vial and allowed to oviposit for 4 days, and then moved again to a new vial and allowed to oviposit for 2 days. Vials were maintained for at least 12 days until all offspring had eclosed. Prior to counting offspring, vials were kept in a freezer. Fertility was defined as the sum of all eclosed offspring produced per female over a maximum timespan of 8 days after a single mating.

Lifespan

Line level measurements

Parsch and Pasyukova Labs

Lifespan (i.e., the duration of adult lifespan, in days) was measured for at least 15 isofemale lines per population, with 5-8 replicate vials per line. Each vial contained 10 flies per sex (20 flies in total), collected within 24 hours of eclosion and kept on 5 mL of food medium. Replicates were assayed in two batches (blocks), with at least one replicate for each line in each block. Food vials were changed and mortality was recorded every Monday, Wednesday, and Friday in the JP lab and daily in the EP lab (except for weekends). Flies that had escaped or which had died from non-natural causes during the course of the experiment were marked as censored for subsequent analysis.

Population level measurements

Flatt Lab

In contrast to the JP and EP labs, the TF lab assayed adult lifespan at the population level, not at the level of isofemale lines nested in populations. 24-hour cohorts of adult flies were kept in 1L demography cages (see Tatar et al., 2001 for details of cage design), with 10 replicate cages per population (each cage contained 5 flies per line and sex for each of the populations; see https://github.com/esradm/DrosEU_PhenotypingWG/tree/main/SummaryTables for a detailed list of lines). Age at death was scored when changing vials (containing 5 mL of food) on the cages, at first every second day for the first 3 weeks of the experiment, and thereafter on Mondays, Wednesdays and Fridays, until all individuals in the experiment had died. Flies that had escaped or which had died from non-natural causes during the course of the assay were marked as censored for subsequent analysis.

Cold-shock mortality

Gonzalez, Kozeretska and Vieira Labs

Measurements were performed in 32, 9 and 2 batches (blocks) for the JV, JG and IK labs, respectively. For each isofemale line and sex, 1-8 replicates consisting of 6-25 flies collected around the peak of eclosion time were placed in vials with fresh food at least one day before the experiment. Five- to seven-day-old flies were then transferred into empty vials immersed in an ice-water bath in a box (polystyrene or styrofoam) and stored at 4°C. After 18 hours of cold-shock, vials were moved to a room at 25°C, and the number of dead flies was scored at a single time point 24 hours later. Cold-shock mortality was estimated for each vial as the proportion of dead flies, i.e. the number of dead flies divided by the number of assayed flies.

Chill-coma recovery time

Vieira and Mensch Labs

Measurements of chill-coma recovery time (CCRT) were performed in 27 and 11 batches (blocks) for the JV and JM labs, respectively. Six 7-day-old flies per isofemale line and sex were placed in an empty vial (one vial per sex per line) immersed in an ice-water bath in a polystyrene box placed in a 4°C room. Six hours later, flies were removed from the vials and placed into individual wells of 24-well plates while being kept on ice. A timer was started once plates were moved from the ice to a bench in a room at 25°C. The recovery of each fly was monitored for a maximum duration of 60 minutes. Flies that were able to stand on their legs were considered recovered and CCRT (in seconds) was recorded. Flies that did not recover within 60 minutes of the recovery period were marked as censored for subsequent analysis.

Heat-shock mortality

Parsch and Vieira Labs

Heat-shock mortality was measured for at least 15 isofemale lines per population. Single sex groups of fifteen 5- to 7-day-old flies were placed into empty vials inside a incubator set at 37ºC and the number of dead flies was scored for 7 hours every 30 minutes in the JV lab, or for 8 - 8.5 hours approximately every 45 minutes (for a total of 10-11 observations) in the JP lab. For each isofemale line, 5 replicates were assayed per sex. Replicates were measured in 32 batches (blocks) for the JV lab, and in 9 batches (blocks) during a single day for the JP lab. Heat-shock mortality was estimated for each vial as the proportion of dead flies, i.e. the number of dead flies divided by the number of assayed flies.

Reproductive diapause

Bergland, Flatt and Schlötterer Labs

To induce adult reproductive diapause (or dormancy), 2-hour-old virgin females (on average 15 flies per isofemale line per population) were exposed to standard diapause-inducing conditions (Saunders et al., 1989), i.e., 12°C and 10:14 hours light:dark, during 3 weeks. Flies were transferred to new vials once per week. After 3 weeks under diapause conditions, flies were kept at -80ºC until dissection. Both ovaries were examined; an individual was classified as diapausing if all oocytes < stage 10 and if no mature eggs were present; individuals were classified as non-diapausing if at least one oocyte was > stage 10 or if mature eggs were present.

Locomotor activity

Tauber Lab

Locomotor activity was measured on 1-13 males (on average 4) for each line using DAM2 Drosophila monitors (Trikinetics Inc., Waltham, MA, USA) in 2 batches (blocks). Single 1- to 3-day-old flies were placed into vials (10 cm x 0.5 cm) filled with 2 cm sugar/agar medium. Monitors were placed in light (LED) chambers, in an incubator at 24°C with ~30% humidity. Flies were entrained to a light-dark cycle (LD 12:12) for 5 days and then allowed to free-run for 10 days in constant darkness (DD).

Circadian eclosion timing

Wegener Lab

Eclosion rhythmicity was measured at the population level using outcrossed isofemale lines for each of the 9 populations (Table XX for IDs of lines used for outcrossing). For each experiment, similar numbers of offspring for all isofemale lines of a given population were pooled and interbred. Flies were raised at a light-dark cycle of 14:10 hours (LD 14:10) in a climate chamber at either 18°C or 29°C, depending on the experiment. Age-mixed puparia of the resulting F1 generation (spanning an age difference of around 5 days at 18°C, or 4 days at 29°C) were collected and glued to a circular perspex disc using fungicide-free methyl cellulose glue (Tapetenkleister Nr. 389, Auro, Germany; diluted 1:30 in water). Eclosion was monitored for one week under LD14:10 or under constant darkness (DD) at either 18°C or 29°C using Drosophila eclosion monitors (Trikinetics Inc., Waltham, MA, USA).

Pigmentation

Abbott, Gibert and Schmidt Labs

For each line, ten 13- to 15-day-old females, either alive or preserved in 95% ethanol, were air dried and placed on their left side, with photographs taken using a dissecting microscope (PG lab: Axio Imager Z1, Zeiss; JA lab: Nikon SMZ1270, PS lab: Olympus DP73 digital camera mounted on a Leica MZ9.5 stereo microscope). Images were analyzed with ImageJ 1.46r, using the “Area Fraction” measurement tool. “Area Fraction” measures the percentage of pixels in a selected area highlighted in red using the “Threshold” tool, yielding an estimate of the percentage of dark pigmentation on the three terminal tergites of the abdomen (tergites 4, 5 and 6). The same tergites were scored by the PS lab using the procedure described in David et al. (1990); pigmentation scores were multiplied by 10 in order for them to be converted into pigmentation percentages similar to those measured by the JA and PG labs.

Starvation resistance

Gonzalez, Onder and Pasyukova Labs

Measurements were performed in 5 batches (blocks) with at least 1 replicate per line in each batch. For each isofemale line and sex, 10 replicate vials, each with ten 3- to 7-day-old flies were assayed. Flies were kept in vials with 5 mL of 2% agar for the duration of the assay. Age at death was scored every 8 hours and starvation resistance was estimated as the number of hours from the start of the experiment until death.

Parasitoid resistance

Hrcek Lab

Parasitoid resistance was measured using 5-8 isofemale lines per population in 3 batches (blocks). Groups of 3- to 5-day-old adults (at least 50 pairs) were allowed to lay eggs on an agar plate overnight for 10 hours. Yeast paste was provided to stimulate egg laying. Immediately thereafter, 80 eggs were collected and placed into a vial with 10 mL of medium. Eight to ten vials were prepared per line, depending on the number of eggs collected. Two days after eggs collection, a single female parasitoid Leptopilina boulardi, aged between 5- to 7-days old, was introduced for 24 hours into each vial (4 to 6 vials) belonging to the parasitization treatment group. Four vials were left non-parasitized as a control. Parasitoid resistance was measured as the proportion of adult flies emerging from the parasitized vials divided by the mean emergence of adults in the non-parasitized control vials.

Wolbachia infections

Isofemale lines were screened for the presence of Wolbachia using three complementary approaches: (1) regular PCR-based screening using protocols described in Miller et al. (1988) and Faria & Sucena (2017) (ES lab); (2) this first screen was independently repeated using regular PCR following Strunov et al. (2022) (MK lab); and (3) quantitative real-time PCR (qPCR)-based screening following a modified protocol from Sambrook et al. (1989) (EP lab). For approach (1), whole-genome DNA was extracted from 3-5 flies per line following the protocol of Miller et al. (1988). A multiplex PCR reaction was carried out in 96-well plates using 1μL of diluted DNA with the GoTaq G2 Flexi DNA polymerase (Promega) in a 10μL total reaction volume per well using the following diagnostic primer pairs: Slif (Fwd: 5’ GTTAGCGCCTATTAGCACAT 3’; Rev: 5’ CGGGACAACTCAGTCTGTAA 3’); wsp (81Fwd: 5’ TGGTCCAATAAGTGATGAAGAAAC 3’; 691Rev: 5’ AAAAATTAAACGCTACTCCA 3’). The latter pair of primers serves to diagnose the presence or absence of Wolbachia (Zhou et al., 1998). The following PCR protocol was used: 95 ºC for 10 min; 30 cycles at 95 ºC for 30s, 60 ºC for 1 min, 72 ºC for 1 min; and a final extension step at 72 ºC for 10 min. PCR products were visualized by gel electrophoresis (1.5% agarose in TAE supplied with 1% RedSafe). For approach (2), we performed PCRs using VNTR-141 primers (Riegler et al., 2012) following the PCR conditions described in Strunov et al. (2022). This allowed us to distinguish between the two most common Wolbachia variants, wMel and wMelCS, based on a diagnostic length polymorphism in the VNTR region. For qPCR, whole-genomic DNA was extracted from pools of 20 flies per line following Sambrook et al. (1989). For approach (3), qPCR was performed on a MiniOpticon real-time PCR system (Bio-Rad). The total reaction volume was 20 µl per well, with 1 µl of diluted DNA, HotStart Taq (Sibenzyme), and SYBR Green I and W-Spec primers (Fwd: 5’ CATACCTATTCGAAGGGATAG 3’; Rev 5’ AGCTTCGAGTGAAACCAATTC 3’; Werren & Windsor, 2000). The W-Spec primer set is thought to give a stronger and more specific signal than a number of other primers (Simoes et al., 2011). The following PCR protocol was used: 94 ºC for 2.5 min; 50 cycles at 94 ºC for 20 s, 64 ºC for 20 s, 72 ºC for 30 s. The specificity of the PCR products detected was determined by melting curve analysis.

1.6 -Statistical analyses

Linear models

Because of (for practical reasons inevitable) differences in data collection and structure among labs, it was typically not feasible for us to fit a single global trait-specific model that would incorporate all data from several labs that had measured the same phenotypic trait. Instead, we had to fit lab-specific models. In total, we ran 97 individual linear models (see XX for the separate analyses of circular data, i.e., circadian eclosion timing and locomotor activity). Modeling was performed in R (v.4.1.1) with the lmer function from the afex package (v.1.0-1). In all models, Population (i.e., population ID) was included as a fixed factor and, whenever applicable, Line (isofemale line ID), Replicate (e.g., replicate vial) and Batch (block) were included as random factors with appropriate nesting as in the following example:
Trait ~ Population + (1|Line:Population) + (1|Batch) + (1|Replicate:Line:Population). To keep the number of factors (and interactions) in the models small, and because we were not specifically interested in quantifying sexual dimorphism, the sexes were analyzed separately. Proportional data (e.g., for traits such as viability, pigmentation and cold-shock mortality) were arcsine-square-root transformed prior to analysis. Diapause and parasitoid resistance data were analyzed with binomial generalized linear mixed-effect models (glmer function from the afex package version 1.0-1) with Population as a fixed effect and Line nested within Population as a random effect. The number of flies scored and the number of flies emerging in control vials served as weights for diapause and parasitoid resistance, respectively. In two cases (locomotor activity ‒ absolute phase measured by the ET lab and viability measured by the PS lab) we used linear (fixed-effects) models instead of linear mixed-effect models because of singularity issues or lack of line replication. For the 95 linear mixed-effect models marginal R2 values (proportion of variance explained by Population) were extracted using the r2_nakagawa function from the performance package (v.0.10.2) (Nakagawa & Schielzeth, 2013). For the two linear models mentioned above we extracted R2 values from the model summaries. Population estimates and associated standard errors were extracted from all 97 model outputs by using the emmeans function from the emmeans package (v.1.7.1-1). Line random coefficients and associated standard errors were derived from Line random effects extracted from model outputs using the extract_random_effects function from the mixedup package (v.0.3.9) (Clark, 2022). Population coefficients were subsequently added to Line random effects to obtain Line random coefficients. Standard errors for the Line random coefficients were calculated as follows: sqrt(Population_SE^2 + Line_random_effect_SE^2)).

We note that in some cases (e.g., for CCRT) Q-Q plots indicated departures from normality. In spite of this, we did not do anything special to deal with such cases as ANOVA/linear models are known to be extremely robust to departures from non-normality. Blanca et al. (2017), for example, conducted a comprehensive, systematic Monte-Carlo simulation study (using a wide variety of common, non-normal distributions) and found that, under a very broad range of conditions (actually in 100% of all simulated cases), non-normality did not present any problems, with F-tests being completely robust in terms of Type 1 error. Thus, we trust/assume that our analyses are likely robust despite non-normality, however, we can of course not prove that non-normality might not have inflated the Type 1 error rate in these cases.

Analyses of survival and mortality data

In addition to using linear mixed-effects models for analyzing time-to-event data (i.e., for lifespan, cold-shock mortality, chill-coma recovery time, heat-shock mortality, starvation resistance), we also used mixed-effects Cox models (i.e., Cox regression; proportional hazard analysis) implemented in the coxme package in R to analyze these traits. ### Analyses of locomotor activity Locomotor activity data were processed into 30 minute bins, and four variables were analyzed. These included (1) circadian period and (2) acrophase (phase in DD) and were analyzed with the MESA algorithm (https://biodare2.ed.ac.uk/) (Zielinski et al., 2014). The other two variables, (3) activity level and (4) nocturnal/diurnal ratio, were analyzed with a R script (Pegoraro et al., 2022). The level of activity represents the daily number of 30-minute bins in which the animal moved at least once, averaged over 5 days. The ND ratio was calculated as the total number of 30-minute bins in which the animal was active during the 12 dark hours, divided by the number of activity bins during the 12 light hours (over 5 days). Except for the acrophase, data analysis was performed in R (R Development Core Team 2013). Except for the acrophase, data were analyzed with mixed-model ANOVA using the lme function in the nlme package in R (Pinheiro et al., 2023). Line was treated as a random effect. Analysis of acrophase was carried out using the Oriana software for circular statistics (Kovach Computing Services, Pentraeth, Isle of Anglesey, UK).

Analyses of circadian eclosion timing and eclosion rhythmicity

Eclosion events were processed into 1 hour bins. Phase and period were calculated using maximum entropy spectrum analysis (MESA); rhythmicity was assessed by using the JTK_CYCLE model, followed by Benjamini-Hochberg correction, as well as using Lomb-Scargle periodogram analysis, implemented in BioDare2 (https://biodare2.ed.ac.uk) (Zielinski et al., 2014).

Effects of Wolbachia infection status on phenotypic variation

Based on the Wolbachia infection status of the isofemale lines as determined by the methods described above, we tested for phenotypic effects of Wolbachia presence in all investigated populations. To do so, we focused on isofemale lines that were unambiguously identified as being Wolbachia-infected (wol+) or uninfected (wol-) with both PCR approaches described above. The populations from Finland and Russia were excluded as all lines from these locations were Wolbachia-infected. Moreover, we only included lines that were assayed in all labs investigating a given phenotype. Finally, we only considered populations with at least three isofemale lines of each infection type (wol+, wol-). For each trait, we fitted linear mixed-effects models using the lme4 package in R (Bates et al., 2015) to test for the fixed effects of Population and of Wolbachia infection status on phenotypic variation. Whenever possible, we also included the factors Sex and Protein to carbohydrate ratio (P:C-ratio) of the laboratory fly food. To account for the fact that the latter factor is confounded by lab identity, we additionally included the random factor Lab in our models. Moreover, we fitted the random factors Line (nested within population) and Experimental batch (nested within lab) to account for biological and technical replication. To test for significance of the fixed factors and all possible interactions, we employed Type-III analysis of deviance using the R package car.

Heritability estimates; genotypic, environmental and phenotypic variances

For each trait and lab, we estimated an upper limit to the broad-sense heritability (H2, i.e., the total genotypic variance divided by the phenotypic variance), the so-called “isofemale heritability” (i.e., the intraclass correlation), by assuming that isofemale lines faithfully represent distinct genotypes (see Parsons, 1983; Hoffmann & Parsons, 1988; Falconer & Mackay, 1996; Lynch & Walsh 1998; David et al., 2005). To do so, we estimated phenotypic and genotypic (as well as environmental) variances using linear mixed-effects models with the lmer function in the lme4 package in R (Bates et al., 2015); the models had the following form: Y = L + e, where Y is the phenotypic trait, L represents the random effect of genotype (i.e., the Line IDs), and e is the error. We stress that these estimates provide an upper bound of H2 and will be inflated by the presence of any common environmental effects (Lynch & Walsh 1998).

Correlations between trait values measured by different labs

To assess the extent of reproducibility (repeatability) of trait estimates across labs, we estimated pairwise Pearson’s correlations between trait values estimated by different labs that had measured the same trait. Input trait values for these analyses consisted of Line random coefficients extracted from linear mixed-effect models (see above); correlations were calculated only for traits that had been measured in more than two labs.

Meta-analyses of population differences

For each trait and sex, we used meta-analysis to compute Population summary effects. Each lab/assay in which the effect of Population on a trait was assessed with a linear mixed-effect model was considered to represent a separate “study”. Because Population has 9 levels, we carried out subgroup meta-analysis by considering each population as a subgroup, thus enabling us to test for phenotypic differences between populations which were unlikely to be caused by differences in environmental/assay conditions among labs. Input data for this analysis consisted of estimates and associated standard errors for the factor Population obtained from trait- and lab-specific linear mixed-effect models. Estimates were used as Population effects and standard errors were used as weights, i.e., to give more or less weight to “studies” (labs) depending on sample sizes and replication levels. Data from the NP lab (wing area, thorax length) were excluded as our threshold for analysis was a minimum of 5 lines per population but as only 3 lines per population had been phenotyped. The exclusion of these data, and the fact that the MR lab could only assay 5 out of 9 populations, left us with data from only a single single lab (IK), preventing us from analyzing male thorax length. Subgroup meta-analyses were performed in R (v.4.1.1) using a random-effects model implemented in the metagen (random = TRUE, method.tau = “REML”) and update.meta (subgroup = Population, tau.common = FALSE) functions of the meta package (version 5.1-1) (Balduzzi et al, 2019). Cochran’s Q was used to assess heterogeneity (i.e., differences in effect sizes) between subgroups (=populations) which were unlikely due to differences in conditions among labs. Resulting p-values were corrected for multiple testing with the Bonferroni procedure (𝛼’ = 𝛼/n = 0.05/26 = 0.0019; n = 26 meta-analyses in total). Population summary effects were extracted from meta-analysis outputs and used as Population compound estimates for downstream analyses. A similar meta-analysis approach (without heterogeneity tests) was employed to generate Line compound estimates using the Line random coefficients and associated standard errors, which were both extracted from the mixed-effect models in which Line could be included as a random effect (see previous section).

Multivariate analyses I: Principal Component Analysis (PCA)

To analyze multivariate phenotypic correlations in more detail, trait values from model estimates were transformed using principal component analysis (PCA) for each sex separately using 134 male lines and 165 female lines (males: 10-17 lines per population, mean = 14.8; females: 14-20 lines per population, mean 18.33). Separate PCAs were carried out for males and females. For males, nine phenotypes were used (chill-coma recovery time, cold-shock mortality, egg-to-adult developmental time, dry weight, heat-shock mortality, lifespan, starvation resistance, thorax length, and wing area (left)). For females a greater number of phenotypic traits were available. To get a more complete picture of multivariate phenotypic correlations, two PCAs with different combinations of phenotypes were carried out: (i) nine phenotypes corresponding to those used in males ‒ see above (i.e., to allow a direct comparison with the male PCA); and (ii) 13 traits including the nine previously mentioned plus three measured only in females (reproductive diapause, fertility, and total pigmentation) and viability (although sex was not specified for viability, it showed strong variation between populations, suggesting that it could be a key phenotype). The male PCA was termed M9 and the two female PCAs were termed F9 and F13. Phenotypic traits were scaled to unit variance and transformed using the PCA function from R package FactoMineR (Lê et al 2008). Traits with greater-than-average contributions to a given PC were considered key correlates of PCs, and phenotypic trait loadings are reported for the first 3 PCs. PCA plots with confidence ellipses were visualized with factoextra (Kassambara and Mundt, 2020). To check for effects of diet, all PCAs were also run with trait values derived from a single lab and chosen such that variation in diet P:C (protein:carbohydrate) ratio was controlled (P:C ratio between 0.066 and 0.283). For males, thorax Length (TL) values were not available at the target P:C ratio, so a complete data PCA and a controlled P:C ratio PCA were run using data from the eight remaining phenotypic traits. Trait loadings for P:C controlled PCAs were compared to trait loadings for the PCAs containing data for all labs using Tucker’s Congruence Coefficients, implemented using the factor.congruence function in the psych package (Revelle, 2022).

Multivariate analyses II: Discriminant Function Analysis (DFA)

Broad differences in population phenotypes were then quantified for each sex, using discriminant function analysis (DFA) and the same 134 male lines and 165 female lines as described above. DFA incorporated the model estimate trait values for 10 phenotypes in males (the nine that were measured plus viability) and 13 in females (the same as the F13 PCA). For each sex, Mahalanobis distance (D2) between populations was calculated, and the observations (134 and 165 lines in males and females, respectively) were reallocated to determine how distinct the populations were. All analyses were run in Genstat v22.1 with Genstat Procedure Library v30.1 loaded using the discriminate function.

Multivariate analyses III: Correlations between environment and phenotype

To analyze effects of environmental factors on phenotypic traits, we downloaded for each population climatic data from the NASA database (https://data.nasa.gov) using geographical coordinates (latitude, longitude) of our sampling locations with help of nasapower R package (Sparks, 2018). Climate data for two different time periods were used: 30 years, in order to test for potential long-term effects on traits, and 30 days prior to the collection (sampling) date in order to test for effects of variation in weather conditions that might have impacted the grandparental or parental populations of field-collected flies. In total, 14 environmental variables were used for each time period (see Results section 2.12) and transformed with PCA in the FactoMineR package (Lê et al., 2008). Environmental PCs with eigenvalues > 1 were then used as predictors in linear models of multivariate phenotypes. Separate models were run for males and females for all phenotypic PCs of interest. We also tested whether correlations between phenotype PCs and environmental PCs were greater than expected by random chance by using a permutation-based approach. In this analysis, we permuted environmental PC values for PC1 and PC2 among population labels and then re-ran correlations between environmental PCs and phenotypic PCs. In addition, in order to assess the clinality of traits, we estimated Pearson’s correlations between traits and latitude, longitude, and altitude; input trait values consisted of Line compound estimates from meta-analyses (see above).

2 - Results

2.1 - Linear Models

Please note that “Plots and Linear Models by Lab” are presented in alphabetical order.

Here is a summary table with p-values for all the analyses below. p-value adjustments are only made on meta analyses (Bonferroni correction for n = 20 traits).

Viability

Contributors

Gibert Lab: Cristina Vieira, Laurence Mouton, Natacha Kremer, Sonia Martinez, Patricia Gibert

Grath Lab: Ingo Müller, Sonja Grath

Hoedjes Lab: Hristina Kostic, Katja Hoedjes

Schmidt Lab: Ozan Kiratli, Yonatan Babore, Liam Forsythe, Paul Schmidt

Stamenkovic-Radak Lab: Marija Savic Veselinovic, Marija Tanaskovic, Aleksandra Patenkovic, Mihailo Jelic, Katarina Eric, Pavle Eric, Slobodan Davidovic, Marina Stamenkovic-Radak

Zwaan Lab: Joost van den Heuvel, Bas Zwaan

str(droseu$via)

## 'data.frame':    2367 obs. of  17 variables:
##  $ Supervisor.PI                    : Factor w/ 6 levels "Gibert","Grath",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Diet                             : Factor w/ 1 level "NS": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Batch                            : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Population                       : Factor w/ 9 levels "AK","GI","KA",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Line                             : Factor w/ 172 levels "AK1","AK10","AK11",..: 1 1 1 2 2 2 4 4 4 9 ...
##  $ ReplicateVialOld                 : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 1 2 3 1 2 3 1 ...
##  $ ReplicateVial                    : Factor w/ 2367 levels "Gibert_1_AK1_1",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ ProportionEggtoAdultSurvival     : num  0.68 0.73 0.63 0.85 0.75 0.8 0.85 0.88 0.7 0.68 ...
##  $ Country                          : Factor w/ 9 levels "Austria","Denmark",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Latitude                         : num  61.1 61.1 61.1 61.1 61.1 61.1 61.1 61.1 61.1 61.1 ...
##  $ Longitude                        : num  23.5 23.5 23.5 23.5 23.5 ...
##  $ Altitude                         : num  88 88 88 88 88 88 88 88 88 88 ...
##  $ Population_Lat                   : Factor w/ 9 levels "YE","RE","GI",..: 9 9 9 9 9 9 9 9 9 9 ...
##  $ Population_Lon                   : Factor w/ 9 levels "RE","GI","KA",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ Population_Alt                   : Factor w/ 9 levels "KA","AK","GI",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ ProportionEggtoAdultSurvival_asin: num  0.97 1.024 0.917 1.173 1.047 ...
##  $ Color                            : chr  "#A00E00" "#A00E00" "#A00E00" "#A00E00" ...

Descriptive Statistics

Descriptive statistics at the line level, with batch information:

Descriptive statistics at the line level, without batch information:

Descriptive statistics at the population level, with batch information:

Descriptive statistics at the population level, without batch information:

Plots and Linear Models by Lab

min_Via <- min(droseu$via$ProportionEggtoAdultSurvival)
max_Via <- max(droseu$via$ProportionEggtoAdultSurvival)

y-axis is scaled by the minimum (0) and maximum (1) values in the full data set.

Gibert Lab

lmers_anova$Via_Gibert_lmer_pop

## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Population 0.80451 0.10056     8 153.66  8.0873 4.389e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmers_sum$Via_Gibert_lmer_pop

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## ProportionEggtoAdultSurvival_asin ~ Population + (1 | Line:Population) +  
##     (1 | Batch)
##    Data: filter(droseu$via, Supervisor.PI == "Gibert")
## 
## REML criterion at convergence: -547.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.1956 -0.5255 -0.0195  0.5213  2.8588 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  Line:Population (Intercept) 0.012426 0.11147 
##  Batch           (Intercept) 0.000125 0.01118 
##  Residual                    0.012435 0.11151 
## Number of obs: 532, groups:  Line:Population, 169; Batch, 3
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    1.06673    0.02957  27.03333  36.074  < 2e-16 ***
## PopulationGI  -0.16344    0.04382 154.67733  -3.730 0.000269 ***
## PopulationKA  -0.04174    0.04065 155.52140  -1.027 0.306077    
## PopulationMA  -0.12876    0.04058 153.42366  -3.173 0.001822 ** 
## PopulationMU  -0.04274    0.04059 154.79040  -1.053 0.293955    
## PopulationRE  -0.12273    0.04311 155.47126  -2.847 0.005017 ** 
## PopulationUM  -0.05047    0.04169 154.54796  -1.211 0.227874    
## PopulationVA  -0.14130    0.04059 154.79040  -3.481 0.000649 ***
## PopulationYE  -0.26929    0.04059 154.79040  -6.634 5.16e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltGI PpltKA PpltMA PpltMU PpltRE PpltUM PpltVA
## PopulatinGI -0.631                                                 
## PopulatinKA -0.681  0.463                                          
## PopulatinMA -0.684  0.463  0.500                                   
## PopulatinMU -0.682  0.463  0.500  0.500                            
## PopulatinRE -0.642  0.437  0.470  0.470  0.471                     
## PopulatinUM -0.664  0.451  0.487  0.487  0.487  0.458              
## PopulatinVA -0.682  0.463  0.500  0.500  0.501  0.471  0.487       
## PopulatinYE -0.682  0.463  0.500  0.500  0.501  0.471  0.487  0.501

Grath Lab

lmers_anova$Via_Grath_lmer_pop

## Type III Analysis of Variance Table with Satterthwaite's method
##              Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
## Population 0.077956 0.038978     2 27.308  1.9446 0.1624

lmers_sum$Via_Grath_lmer_pop

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ProportionEggtoAdultSurvival_asin ~ Population + (1 | Line:Population)
##    Data: filter(droseu$via, Supervisor.PI == "Grath")
## 
## REML criterion at convergence: -123.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1670 -0.5714  0.0117  0.4894  2.9754 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  Line:Population (Intercept) 0.004157 0.06447 
##  Residual                    0.020045 0.14158 
## Number of obs: 147, groups:  Line:Population, 30
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   0.98774    0.02906 28.44269  33.995   <2e-16 ***
## PopulationMU -0.07674    0.04075 27.59627  -1.883   0.0703 .  
## PopulationRE -0.05958    0.04075 27.59627  -1.462   0.1551    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltMU
## PopulatinMU -0.713       
## PopulatinRE -0.713  0.508

Hoedjes Lab

lmers_anova$Via_Hoedjes_lmer_pop

## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq  Mean Sq NumDF DenDF F value    Pr(>F)    
## Population 0.51582 0.064478     8   158  6.2599 4.985e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmers_sum$Via_Hoedjes_lmer_pop

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ProportionEggtoAdultSurvival_asin ~ Population + (1 | Line:Population)
##    Data: filter(droseu$via, Supervisor.PI == "Hoedjes")
## 
## REML criterion at convergence: -549.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3516 -0.5255 -0.0429  0.4980  4.2544 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  Line:Population (Intercept) 0.01545  0.1243  
##  Residual                    0.01030  0.1015  
## Number of obs: 501, groups:  Line:Population, 167
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    1.09520    0.03073 157.99999  35.644  < 2e-16 ***
## PopulationGI  -0.20137    0.04693 157.99999  -4.290 3.10e-05 ***
## PopulationKA  -0.06612    0.04345 157.99999  -1.522 0.130106    
## PopulationMA  -0.08792    0.04345 157.99999  -2.023 0.044720 *  
## PopulationMU  -0.04924    0.04345 157.99999  -1.133 0.258863    
## PopulationRE  -0.17978    0.04693 157.99999  -3.831 0.000184 ***
## PopulationUM  -0.10841    0.04533 157.99999  -2.392 0.017952 *  
## PopulationVA  -0.11483    0.04345 157.99999  -2.643 0.009052 ** 
## PopulationYE  -0.24862    0.04345 157.99999  -5.722 5.15e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltGI PpltKA PpltMA PpltMU PpltRE PpltUM PpltVA
## PopulatinGI -0.655                                                 
## PopulatinKA -0.707  0.463                                          
## PopulatinMA -0.707  0.463  0.500                                   
## PopulatinMU -0.707  0.463  0.500  0.500                            
## PopulatinRE -0.655  0.429  0.463  0.463  0.463                     
## PopulatinUM -0.678  0.444  0.479  0.479  0.479  0.444              
## PopulatinVA -0.707  0.463  0.500  0.500  0.500  0.463  0.479       
## PopulatinYE -0.707  0.463  0.500  0.500  0.500  0.463  0.479  0.500

Schmidt Lab

lmers_anova$Via_Schmidt_lm_pop

## Analysis of Variance Table
## 
## Response: ProportionEggtoAdultSurvival_asin
##             Df  Sum Sq  Mean Sq F value   Pr(>F)   
## Population   8  1.7653 0.220666  2.6999 0.008308 **
## Residuals  153 12.5050 0.081732                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmers_sum$Via_Schmidt_lm_pop

## 
## Call:
## lm(formula = ProportionEggtoAdultSurvival_asin ~ Population, 
##     data = filter(droseu$via, Supervisor.PI == "Schmidt"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73241 -0.21504  0.01621  0.15661  0.83839 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.92028    0.06393  14.396   <2e-16 ***
## PopulationGI -0.18787    0.09765  -1.924   0.0562 .  
## PopulationKA  0.21332    0.09041   2.360   0.0196 *  
## PopulationMA  0.01281    0.09288   0.138   0.8905    
## PopulationMU -0.03180    0.09041  -0.352   0.7255    
## PopulationRE  0.02327    0.09765   0.238   0.8120    
## PopulationUM  0.10329    0.09962   1.037   0.3015    
## PopulationVA -0.03068    0.09041  -0.339   0.7348    
## PopulationYE -0.08020    0.09041  -0.887   0.3764    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2859 on 153 degrees of freedom
## Multiple R-squared:  0.1237, Adjusted R-squared:  0.07789 
## F-statistic:   2.7 on 8 and 153 DF,  p-value: 0.008308

Stamenkovic-Radak Lab

lmers_anova$Via_StamenkovicRadak_lmer_pop

## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq  Mean Sq NumDF  DenDF F value    Pr(>F)    
## Population 0.50121 0.062652     8 155.29  5.2001 9.104e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmers_sum$Via_StamenkovicRadak_lmer_pop

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## ProportionEggtoAdultSurvival_asin ~ Population + (1 | Line:Population) +  
##     (1 | Batch)
##    Data: filter(droseu$via, Supervisor.PI == "StamenkovicRadak")
## 
## REML criterion at convergence: -485.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0263 -0.5091 -0.0167  0.5112  3.1782 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  Line:Population (Intercept) 0.015859 0.1259  
##  Batch           (Intercept) 0.001318 0.0363  
##  Residual                    0.012048 0.1098  
## Number of obs: 501, groups:  Line:Population, 167; Batch, 4
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    0.98823    0.03644  24.11360  27.117  < 2e-16 ***
## PopulationGI  -0.09533    0.05031 155.41420  -1.895  0.05998 .  
## PopulationKA  -0.04514    0.04464 155.26320  -1.011  0.31344    
## PopulationMA  -0.14459    0.04460 155.11506  -3.242  0.00145 ** 
## PopulationMU   0.01274    0.04460 155.09434   0.286  0.77548    
## PopulationRE  -0.07059    0.04737 155.37260  -1.490  0.13819    
## PopulationUM  -0.10798    0.04589 155.37865  -2.353  0.01987 *  
## PopulationVA  -0.13926    0.04460 155.11479  -3.122  0.00214 ** 
## PopulationYE  -0.21115    0.04464 155.26320  -4.730    5e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltGI PpltKA PpltMA PpltMU PpltRE PpltUM PpltVA
## PopulatinGI -0.542                                                 
## PopulatinKA -0.614  0.443                                          
## PopulatinMA -0.611  0.442  0.499                                   
## PopulatinMU -0.613  0.442  0.500  0.500                            
## PopulatinRE -0.579  0.415  0.472  0.471  0.472                     
## PopulatinUM -0.598  0.429  0.488  0.486  0.487  0.460              
## PopulatinVA -0.612  0.445  0.500  0.499  0.500  0.470  0.485       
## PopulatinYE -0.614  0.443  0.501  0.499  0.500  0.472  0.488  0.500

Zwaan Lab

lmers_anova$Via_Zwaan_lmer_pop

## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Population 1.6727 0.20908     8 150.05  6.3093 4.877e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmers_sum$Via_Zwaan_lmer_pop

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ProportionEggtoAdultSurvival_asin ~ Population + (1 | Line:Population)
##    Data: filter(droseu$via, Supervisor.PI == "Zwaan")
## 
## REML criterion at convergence: -124.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5670 -0.4679  0.0236  0.4733  2.8902 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  Line:Population (Intercept) 0.01407  0.1186  
##  Residual                    0.03314  0.1820  
## Number of obs: 524, groups:  Line:Population, 169
## 
## Fixed effects:
##                Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    1.072023   0.035162 146.338153  30.488  < 2e-16 ***
## PopulationGI  -0.134116   0.053744 145.530228  -2.495   0.0137 *  
## PopulationKA   0.006813   0.049907 147.722363   0.137   0.8916    
## PopulationMA  -0.045478   0.049950 148.488267  -0.910   0.3640    
## PopulationMU  -0.053131   0.050371 151.671992  -1.055   0.2932    
## PopulationRE  -0.122316   0.052526 144.039234  -2.329   0.0213 *  
## PopulationUM   0.033358   0.051199 147.170714   0.652   0.5157    
## PopulationVA  -0.019751   0.050311 151.428503  -0.393   0.6952    
## PopulationYE  -0.258184   0.050259 151.942236  -5.137 8.44e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltGI PpltKA PpltMA PpltMU PpltRE PpltUM PpltVA
## PopulatinGI -0.654                                                 
## PopulatinKA -0.705  0.461                                          
## PopulatinMA -0.704  0.461  0.496                                   
## PopulatinMU -0.698  0.457  0.492  0.491                            
## PopulatinRE -0.669  0.438  0.472  0.471  0.467                     
## PopulatinUM -0.687  0.449  0.484  0.483  0.479  0.460              
## PopulatinVA -0.699  0.457  0.492  0.492  0.488  0.468  0.480       
## PopulatinYE -0.700  0.458  0.493  0.492  0.488  0.468  0.480  0.489

Linear Models / Residuals

Gibert Lab

Grath Lab

Hoedjes Lab

Schmidt Lab

Stamenkovic-Radak Lab

Zwaan Lab

Developmental time

Egg-to-pupa developmental time

Contributors

Schmidt Lab: Paul Schmidt

Data Reformatting

str(droseu$dtp)

## 'data.frame':    3391 obs. of  17 variables:
##  $ Supervisor.PI   : Factor w/ 1 level "Schmidt": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Diet            : Factor w/ 1 level "NS": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Batch           : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Population      : Factor w/ 9 levels "AK","GI","KA",..: 8 8 8 8 8 8 8 8 8 8 ...
##  $ Line            : Factor w/ 161 levels "AK1","AK10","AK11",..: 131 131 131 131 131 131 131 131 131 131 ...
##  $ ReplicateVialOld: Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ReplicateVial   : Factor w/ 161 levels "Schmidt_1_AK1_1",..: 131 131 131 131 131 131 131 131 131 131 ...
##  $ Individual      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ DT_EggPupa      : num  120 120 120 120 136 136 136 136 136 136 ...
##  $ Country         : Factor w/ 9 levels "Austria","Denmark",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ Latitude        : num  58 58 58 58 58 ...
##  $ Longitude       : num  33.2 33.2 33.2 33.2 33.2 ...
##  $ Altitude        : num  217 217 217 217 217 217 217 217 217 217 ...
##  $ Population_Lat  : Factor w/ 9 levels "YE","RE","GI",..: 8 8 8 8 8 8 8 8 8 8 ...
##  $ Population_Lon  : Factor w/ 9 levels "RE","GI","KA",..: 9 9 9 9 9 9 9 9 9 9 ...
##  $ Population_Alt  : Factor w/ 9 levels "KA","AK","GI",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ Color           : chr  "#095888" "#095888" "#095888" "#095888" ...

Descriptive Statistics

# Note that the trait has been phenotyped only in Schmidt lab and in one batch.

Descriptive statistics at the line level:

Descriptive statistics at the population level:

Plots and Linear Models by Lab

min_DT_P <- min(droseu$dtp$DT_EggPupa)
max_DT_P <- max(droseu$dtp$DT_EggPupa)

y-axis is scaled by the minimum (96) and maximum (192) values in the full data set.

Schmidt Lab

lmers_anova$DT_P_Schmidt_lmer_pop

## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Population 2515.1  314.38     8 147.95  2.6412 0.009804 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmers_sum$DT_P_Schmidt_lmer_pop

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DT_EggPupa ~ Population + (1 | Line:Population)
##    Data: filter(droseu$dtp, Supervisor.PI == "Schmidt")
## 
## REML criterion at convergence: 26303.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9734 -0.6092 -0.0831  0.3912  4.7796 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  Line:Population (Intercept) 158.4    12.59   
##  Residual                    119.0    10.91   
## Number of obs: 3391, groups:  Line:Population, 161
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  127.5020     2.8668 143.5783  44.475  < 2e-16 ***
## PopulationGI  12.9203     4.7140 157.7547   2.741  0.00684 ** 
## PopulationKA  -0.2098     4.0656 144.7480  -0.052  0.95892    
## PopulationMA   4.5016     4.1179 144.9709   1.093  0.27613    
## PopulationMU   2.2279     4.0527 143.3696   0.550  0.58337    
## PopulationRE   7.7677     4.4068 146.8919   1.763  0.08004 .  
## PopulationUM   3.6588     4.3705 142.4818   0.837  0.40391    
## PopulationVA  13.2633     4.1248 145.8807   3.215  0.00160 ** 
## PopulationYE   7.6982     4.0711 145.8355   1.891  0.06062 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) PpltGI PpltKA PpltMA PpltMU PpltRE PpltUM PpltVA
## PopulatinGI -0.608                                                 
## PopulatinKA -0.705  0.429                                          
## PopulatinMA -0.696  0.423  0.491                                   
## PopulatinMU -0.707  0.430  0.499  0.492                            
## PopulatinRE -0.651  0.396  0.459  0.453  0.460                     
## PopulatinUM -0.656  0.399  0.463  0.457  0.464  0.427              
## PopulatinVA -0.695  0.423  0.490  0.484  0.492  0.452  0.456       
## PopulatinYE -0.704  0.428  0.497  0.490  0.498  0.458  0.462  0.489

Egg-to-adult developmental time

Contributors

Gibert Lab: Cristina Vieira, Laurence Mouton, Natacha Kremer, Sonia Martinez, Patricia Gibert

Grath Lab: Ingo Müller, Sonja Grath

Hoedjes Lab: Hristina Kostic, Katja Hoedjes

Schmidt Lab: Ozan Kiratli, Yonatan Babore, Liam Forsythe, Paul Schmidt

Stamenkovic-Radak Lab: Marija Savic Veselinovic, Marija Tanaskovic, Aleksandra Patenkovic, Mihailo Jelic, Katarina Eric, Pavle Eric, Slobodan Davidovic, Marina Stamenkovic-Radak

Zwaan Lab: Joost van den Heuvel, Bas Zwaan

Data Reformatting

str(droseu$dta)

## 'data.frame':    57609 obs. of  18 variables:
##  $ Supervisor.PI   : Factor w/ 6 levels "Gibert","Grath",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Diet            : Factor w/ 1 level "NS": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Batch           : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Population      : Factor w/ 9 levels "AK","GI","KA",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Line            : Factor w/ 171 levels "AK1","AK10","AK11",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Sex             : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ReplicateVialOld: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ReplicateVial   : Factor w/ 2300 levels "Gibert_1_AK1_1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Individual      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ DT_EggAdult     : num  202 202 202 202 202 202 202 202 202 202 ...
##  $ Country         : Factor w/ 9 levels "Austria","Denmark",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Latitude        : num  61.1 61.1 61.1 61.1 61.1 61.1 61.1 61.1 61.1 61.1 ...
##  $ Longitude       : num  23.5 23.5 23.5 23.5 23.5 ...
##  $ Altitude        : num  88 88 88 88 88 88 88 88 88 88 ...
##  $ Population_Lat  : Factor w/ 9 levels "YE","RE","GI",..: 9 9 9 9 9 9 9 9 9 9 ...
##  $ Population_Lon  : Factor w/ 9 levels "RE","GI","KA",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ Population_Alt  : Factor w/ 9 levels "KA","AK","GI",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Color           : chr  "#A00E00" "#A00E00" "#A00E00" "#A00E00" ...

Descriptive Statistics

Descriptive statistics at the line level, with batch information:

Descriptive statistics at the line level, without batch information:

Descriptive statistics at the population level, with batch information:

Descriptive statistics at the population level, without batch information:

Plots and Linear Models by Lab

min_DT_A <- min(droseu$dta$DT_EggAdult)
max_DT_A <- max(droseu$dta$DT_EggAdult)

y-axis is scaled by the minimum (150) and maximum (394) values in the full data set.