A reigning paradigm in biology is that short-term evolution can be predicted from measures of genetic variation within populations, but that the accuracy of such predictions should decay with time. Here, we show that intrinsic developmental variability and standing genetic variation in wing shape of the two flies, Drosophila melanogaster and Sepsis punctum, are tightly aligned and predict deep divergence in the dipteran phylogeny, spanning >900 taxa and 185 My of evolution. This finding is hard to reconcile with constraint hypotheses invoking a lack of genetic variation as the reason for slow-evolving wing traits unless most of the observed variability is associated with deleterious side effects and effectively unusable for evolution. However, phenotyping of 71 genetic lines of S. punctum revealed no association between variation in wing shape and fitness correlates unrelated to flight, lending no credence to this hypothesis. We also find no evidence for genetic constraints on the pace of wing shape evolution along individual branches of the phylogeny. Instead, correlational selection related to allometric scaling, simultaneously shaping both developmental bias and deep divergence in fly wings, emerges as the most plausible explanation for the observed patterns. This suggests that past forces of selection have shaped the developmental architecture of the dipteran wing such that its long-term evolution can be predicted from its intrinsic variability. These findings challenge our understanding of the fundamental processes governing the emergence of phenotypic variation and its evolution.

Quantifying wing shape and divergence across fly taxa 

We focus on the evolution of wing shape within the Eremoneura, a clade within the Brachycera characterized by the presence of three larval instars. This clade is about 185 million years old and includes the dance and long-legged flies (Empidoidea) as well as the Cyclorrhapha (flies that pupate within the cuticle of the last larval instar (i.e., the puparium)).

To quantify the morphological variation within and between families, we took advantage of illustrations and pictures of fly wings from the taxonomic and systematic literature. An initial data set was sourced from the Manuals of Nearctic Diptera and the Manuals of Afrotropical Diptera . We focused on those families that are represented in the phylogenetic hypothesis generated by Bayless et al.. Additional pictures and illustrations were collected from a wide range of publications.

We only included observations where the location of all 11 two-dimensional landmarks used in Rohner & Berger could be assigned. In total, we collected wing shape data from 827 individuals belonging to 53 families. Using tpsDig2, we manually quantified wing morphology as depicted on the illustrations and images. Additional morphometric data originally collected from images was added for 119 species of drosophilids from Houle et al. 2017 and 36 species for sepsids from Rohner & Berger 2023. The final data set contained 993 observations of 933 species in 530 genera and 68 families. Family affiliation of individual genera was checked using Systema Dipterorum.

The number of observations varied strongly across families (mean = 17.40, median = 12, minimum = 1, maximum = 130). This uneven sampling was caused by i) a varying number of species per family (e.g., Australimyzidae is a monogeneric family containing just 9 described species compared to Tachinidae with 9,626 species), the loss of landmarks in several species (e.g., Sphaeroceridae), and often incomplete illustrations or pictures showing only part of the wing (e.g., Muscidae and Calliphoridae). The landmark coordinates were aligned to the mean configuration of Houle et al. 2017 using Procrustes analysis in MorphoJ.

Quantifying variation in genetic quality to test for deleterious pleiotropy associated with wing shape

To quantify fitness variance across the same isofemale lines of S. punctum as measured for wing shape, we reared all 71 lines (originating from the 7 European populations) in a common garden experiment including 9 temperature treatments ranging from 15-31°C. Five fitness-correlates were measured. F0 containers with fly cultures were seeded with vials of cow dung to attain freshly laid eggs. Each line was seeded with four vials per temperature treatment. For each vial, juvenile development rate was estimated as the inverse of the time (in days) between the date of a laid clutch and the subsequent emergence of F1 adults. Juvenile survival was calculated at the fraction of laid eggs that emerged as adults. Emerged F1 females were measured for their tibia length as an estimate of body size and then paired with a male from the same line and placed in a 50ml vial with access to sugar, water, and cow dung as egg laying substrate. Sugar, water, and dung was replaced every 5 days for the first 15 days. Early reproductive success was estimated as the total number of offspring produced within the first 15 days of adult female life, excluding females that died during this timeframe (likely due to accidental deaths). Females that did not lay any eggs during this period were not included in estimates of early reproductive success. Lifespan was estimated as the time from the start of the experiment until the focal female died. 245 females (17%) did not die during the period of observation and were recorded as censused data. In total, 1,445 females were measured across all lines and temperature treatments. In total, these females produced 173,556 offspring.

To test for heritable variation in early reproductive success, juvenile survival, developmental rate, and body size, we used mixed effects models using restricted maximum likelihood as implemented in ASReml-R . Temperature treatment, population, as well as their interaction were fitted as fixed effect. Line was added as random effect. Note that we did not estimate line by treatment interactions (i.e. G-by-E) as our aim was to capture overall differences in genetic quality among lines. Repeating the analysis excluding the highest (31℃) and lowest (15℃) (i.e. most stressful) temperature treatments led to similar results. Residual variances were allowed to vary across treatments. The significance of the random effect of line was tested using LRTs. Best linear unbiased predictors (BLUPs) were extracted and used for further analysis. For the analysis of adult lifespan, we fitted a censored mixed effects cox model using the coxme package [90] using treatment and population as fixed effects and line as random effect. Likelihood Ratio Tests (LRTs) were used to test whether line effects were significant. We extracted hazard ratios for each line and used its inverse as our estimates for adult longevity.