A cis-regulatory sequence of the selector gene vestigial drives the evolution of wing scaling in Drosophila species

Genetic and environmental factors contribute to morphological diversity across related species (Carroll, 2000). For instance, evolution of morphogenetic traits may be driven by environmental cues followed by selection of specific genetic variations in a population (Uller et al., 2018). The search for genetic variations that account for phenotypic changes across species has attracted the attention of evolutionary biologists for decades, but only recently have genome-editing technologies enabled us to test some predictions directly (Ryan and Farley, 2020). One way to study the genetic contributions to diversity in related species is to investigate how changes in the genome could explain the evolution of certain organs to whole-body scaling relationships (Evans, 2017; Van Belleghem et al., 2020).

To investigate scaling of a specific trait or body part, it is useful to evaluate quantitatively the relationship between the trait and the whole body (or representative trait) within the same individual. When body proportions do not change with body size, simple proportions may be used. Otherwise, scaling relationships should be used to analyze morphological changes that give rise to diversity and phenotypic variation in ecology and evolution (Esquerré et al., 2017; Galicia-Mendoza et al., 2021; Gomes Rodrigues et al., 2018; Shingleton, 2010). Depending on how scaling relationships are studied, allometries can be classified as ontogenic, when scaling relationships are studied during development and growth of an organism (Esquerré et al., 2017; Simons and Frost, 2021); static, when scaling relationships are studied within a population of individuals in a specific stage of development (Shingleton et al., 2008); and evolutionary, when scaling relationships are studied across different species (Tidière et al., 2017).

In Eqn 1, the transformed variables log(x) and log(y) are fitted to a linear transformation using standard linear-fitting procedures, where a is the slope of the line (also known as the allometric coefficient) and log(b) is the intercept (Gayon, 2000; Shingleton, 2010). Using this formalism, allometries have been used to understand how environmental conditions (such as temperature, nutrition and rearing density) influence natural variation and trigger variability in scaling patterns (Shingleton et al., 2009; Weber, 1990), and to investigate the relationships underlying morphological evolution across different species (Pélabon et al., 2014; Shingleton, 2010). However, the genetic mechanisms underlying the evolution of allometric relationships in a group of species are little understood.

One way to understand the genetic basis underlying the evolution and diversity of body traits is to investigate the influence of key genes (e.g. selector genes) on the allometric relationships (Eqn 1). Genetic changes that affect the sensitivity to environmental factors generally affect slopes or allometric coefficients (Emlen et al., 2012; Okada et al., 2019; Tang et al., 2011; Xu et al., 2015), but genetic modifications may also control the intercept of the scaling relationships, such as the wing–body versus haltere–body relationship in Drosophila melanogaster (Crickmore and Mann, 2006). Therefore, changes in the intercept of an organ to the whole-body scaling relationship across species may be attributed to the evolution of specific developmental genes or pathways, such as the regulation of selector genes (Fig. 1). In particular, the evolution of relative organ size may be studied by investigating genetic changes affecting regulatory networks that control the expression pattern of an organ-specific selector gene (Fig. 1A), the properties of the protein that is produced by the selector gene (e.g. stability, degradation or its ability to bind DNA or other proteins; Fig. 1B), or by modifications in cis-regulatory elements (promoters, enhancer or silencers) that drive its expression (Gracia-Latorre et al., 2022; Rebeiz and Tsiantis, 2017; Ryan and Farley, 2020) (Fig. 1C).

The Drosophila wing has been extensively used as a model system to understand the genetic basis of organ growth and morphogenesis (Gou et al., 2020; Neto-Silva et al., 2009; Tripathi and Irvine, 2022). In particular, it is a useful system to investigate quantitatively the variation of body traits under genetic and environmental manipulations. Adult wing size in Drosophila is determined by growth of a precursor organ, the wing imaginal disc, which is patterned and grows during larval development (Beira and Paro, 2016). Growth of the wing imaginal disc is driven by two morphogens: Wingless (Wg) emanating at the dorsal–ventral boundary, and Decantaplegic (Dpp), at the anterior–posterior boundary in the wing imaginal disc (Couso et al., 1995; Neumann and Cohen, 1996; Prasad et al., 2003). These morphogens act together with vestigial (vg), an evolutionarily conserved wing selector gene in holometabolous insects (Abouheif and Wray, 2002; Clark-Hachtel et al., 2013; Nel et al., 2013; Williams et al., 1991) in a regulatory network that triggers vg expression and wing growth through proliferation and cell recruitment (Muñoz-Nava et al., 2020; Parker and Struhl, 2020; Zecca and Struhl, 2007a,b). Through cell recruitment, Vg-expressing cells induce the propagation of the vg pattern via the Quadrant Enhancer (vg^QE) in a complex signal processing known as ‘feed forward signal’, which is fueled by Wg and Dpp signaling and implicates the polarization of the protocadherins Fat and Dachsous, that ultimately results in the nuclear translocation of the transcription factor Yorkie (Yki), where it finally activates transcription of vg (Muñoz-Nava et al., 2020; Neumann and Cohen, 1996; Parker and Struhl, 2020; Zecca and Struhl, 2007a,b, 2010).

Here, we investigated how wing sizes scale relative to tibia sizes (used as a proxy for whole-animal size) in four Drosophila species: D. melanogaster, D. simulans, D. ananassae and D. virilis (Fig. 2A). Although there are some significant deviations in the slopes of the static allometric relationships of these species, we immediately noticed that D. virilis has a dramatically different intercept than the other three species. This difference results from the fact that D. virilis individuals have a much smaller relative wing size compared with the other species. Thus, we asked whether the evolution of this scaling relationship is driven by changes in the cis-regulatory elements within the vg selector gene (Fig. 1C). To address this, we cloned the sequence of the predicted vg^QE in D. virilis and used the CRISPR/Cas9 system to replace these sequences into the D. melanogaster genome. Strikingly, in animals where the vg^QE of D. melanogaster was replaced by the corresponding vg^QE of D. virilis, wings exhibit a reduced relative wing size with respect to what would be expected in D. melanogaster. These results suggest that the vg^QE sequence plays an important role in the determination of size and the establishment of wing-to-body scaling relationships in Drosophila species.

The following Drosophila strains were used: D. melanogaster Meigen 1830 Samarkand strain (RRID:BDSC_4270, Bloomington Drosophila Stock Center); y,w strain of D. melanogaster, provided by Dr Fanis Missirlis (Cinvestav, Mexico); D. simulans Sturtevant 1919 (14021-0251.261, Drosophila Species Stock Center); D. ananassae Doleschall 1858 (14024-0371.00, Drosophila Species Stock Center); D. virilis Sturtevant 1916 (15010-1051.87, Drosophila Species Stock Center); w, ms1096-Gal4 (Bloomington Drosophila Stock Center RRID:BDSC_8860); UAS-vg^RNAi (Vienna Drosophila Resource Center 16896); nub-Gal4 (Bloomington Drosophila Stock Center RRID:BDSC_38418); w; Sco/CyO, RFP (available from Timothy A. Evans); y,w, nos-Cas9 (Bloomington Drosophila Stock Center RRID:BDSC_54591); and D. melanogaster (vg^QEDmel) and D. melanogaster (vg^QEDvir), both generated in this study. Adult flies were crossed and kept at 25°C in vials containing standard Drosophila food.

Adults were separated by sex (data in Figs 2 and 4 correspond to females, but equivalent results were obtained for males in Figs S2 and S4, respectively) using a stereoscopic microscope (Nikon SMZ800) and preserved in 1 ml of 70% ethanol for dehydration for at least 12 h for further analysis. Each specimen was dissected in 15 μl of ethanol 50% to obtain left/right wings, left/right fore legs and thorax, which were mounted in microscope slides using a stereoscopic microscope (Nikon SMZ800).

Wings, fore legs and thorax were photographed in a bright-field binocular microscope (Nikon Eclipse Ci) attached to a camera (ProgRes^® CT5, Jenoptik) using the software ProgRes^® Capture Pro-2.9. Measurements of proximal-distal length (P/D) (see Fig. S2A) were performed using ImageJ software (https://imagej.nih.gov/ij/download.html) (Schneider et al., 2012) after the corresponding calibration for the 4× objective (distance in pixels: 100.501; known distance: 0.1; pixel aspect ratio: 1.0, unit of length: mm). Each point in the allometry graphs corresponds to the average measurements of left/right wings and left/right fore legs of the same animal. We used standardized major axis (SMA) regression (Shingleton, 2019) to obtain the slope and intercepts for each population of flies.

Synteny of the vg gene in D. melanogaster, D. simulans, D. ananassae and D. virilis was evaluated through information from the NCBI (https://www.ncbi.nlm.nih.gov/) and OrthoDB databases (https://www.orthodb.org/). The search was limited to the fourth intron of the gene, and posteriorly paired/local alignments using the FASTA format through the Smith–Waterman algorithm (Smith and Waterman, 1981) in EMBOSS Water (Ubuntu 18.04; https://launchpad.net/ubuntu/bionic/+package/emboss) were developed (D. melanogaster QE versus Drosophila species introns): the length of each sequence was delimited based on the D. melanogaster reference sequence, followed by multiple alignments using Clustal Omega from EMBL-EBI (Madeira et al., 2022). The sequence of QE reported for D. melanogaster strain Samarkand (GenBank ID: FJ513071.1) was used as a template for searching corresponding QE sequences for D. simulans, D. ananassae and D. virilis. Alignment of representative QE regions in Fig. 3 and whole multiple sequence alignment was developed using BioEdit version 7.7.1 (see Fig. S5).

Target sites were designed using the flyCRISPR online tool CRISPR Optimal Target Finder (https://flycrispr.org) (Gratz et al., 2014; Iseli et al., 2007). Cloning was performed into the pCFD4-U6:1_U6:3tandemgRNAs vector Addgene plasmid 49411, which allows tandem expression of two gRNA sequences (Evans, 2017; Port et al., 2014), through PCR amplification using pair of primers 2 and 3 (Table S1) with a 2X Phusion Flash PCR Master Mix (Thermo Fisher Scientific, catalog no. F548S). Gibson Assembly (New England Biolabs, catalog no. E2611) was performed with PCR products and pCFD4 BbsI (New England Biolabs, catalog no. R0539L) digested vector. All cloning products were confirmed by DNA sequencing before injection.

The donor plasmid was assembled using overlap extension PCR (Heckman and Pease, 2007). We first designed a working base vector composed of four PCR fragments: two that were derived from plasmid backbone of pHD-DsRed (Addgene plasmid 51434; https://www.addgene.org/51434/) (primer pair 835, 836), DsRed fluorescent coding region (primer pair 837, 838), and two derived from genomic DNA of D. melanogaster, left homologous arm (LHA2) (primer pair 839, 840) and right homologous arm (RHA1) (primer pair 841, 842). We then assembled the LHA2-DsRed-RHA1 fragments in one round of PCR (primer pair 839, 842) and then combined this product with the plasmid backbone fragment in a second PCR to obtain the whole base vector as a linear fragment: pHD-DsRed-LHA2-DsRed-RHA1 (primer pair 837, 843). The PCR product was transformed directly into competent E. coli and circularized in vivo (Watson and García-Nafría, 2019). To insert the specific QE sequences, we amplified by PCR the sequences from genomic preparations of D. melanogaster (primer pair 844, 845) and D. virilis (primer pair 844, 846) and assembled with the base vector through overlap extension PCR. All the components of the donor vector, including the QE, were sequenced prior to injection.

The vg^QE gRNA plasmid was co-injected with the QE homologous donor plasmids into nos-Cas9 embryos using the service provided by Rainbow Transgenic Flies (Camarillo, CA, USA). Individuals recovered from injected embryos were crossed to w; Sco/CyoRFP flies. CyORFP or Sco flies from this cross (potentially carrying vg^QE HDR alleles) were screened for red fluorescent eyes (indicating genomic integration of DsRed sequences carried on the HDR donor plasmid; see Fig. S3). Fluorescent-eye flies were then crossed back individually to w; Sco/CyORFP to generate balanced stocks. Confirmation of the inserted sequences was tested by PCR using primers 20, 23, 41, 44, 47 and 53 (Table S1, Fig. S3). PCR-confirmed balanced stocks with the modified alleles were also confirmed through genomic DNA sequencing.

Scaling relationship data were analyzed using R version 4.1.2 (https://www.r-project.org/) and RStudio (https://www.rstudio.com/) version 4.1.2. The collected data were analyzed using an SMA model II regression (Shingleton, 2019; Warton et al., 2006; Wilcox et al., 2022 preprint) and computed 95% confidence intervals through the SMATR package in R (Warton et al., 2012; script is available in GitHub). Wing and tibia lengths were analyzed using GraphPad Prism 8.0.1 Software (https://www.graphpad.com/scientific-software/prism/). For pairwise statistical comparisons, Student’s t-tests were used with P<0.05 to define significance.

Although several measures of whole-body size (e.g. mass or pupal case) have been used in Drosophila studies (Mirth and Shingleton, 2012; Morin et al., 1996; Stillwell et al., 2011), we looked for a trait that could serve as a reference to draw allometric relationships with the wings. The length of the tibia in legs has been used as a proxy for body size measures in several dipterans (Krause et al., 2019; Reigada and Godoy, 2005; Rohner et al., 2018; Tran et al., 2020), and is considered a standard for body size in other insects, such as parasitoid wasps (Sagarra et al., 2001; Thorne et al., 2006). Tibias can be much more reliably measured and directly compared with the proximal–distal (PD) length of the adult wing than other whole-body measures such as the thorax (Fig. 2B; see Materials and Methods). In addition, legs are not affected by changes in vg expression (Fig. S1A,B). To further justify the use of tibias as a proxy of body size, we compared the wings, fore tibias and thorax of individuals in each of the four species shown in Fig. 2A (D. melanogaster, D. simulans, D. ananassae and D. virilis; Fig. 2B). Using mean trait measurements in each species (e.g. mean_thorax_Species), we can compute wing, fore tibia and thorax ratios of each species relative to D. melanogaster (e.g. mean_thorax_Species/mean_thorax_Dmel). The thorax ratio (T) reveals that D. simulans is smaller than D. melanogaster (T<1, blue bar of Fig. 2B and Fig. S2A) and D. ananassae is slightly larger (T>1, green bar in Fig. 2B and Fig. S2A), but D. virilis is much larger (T>>1, magenta bar in Fig. 2B and Fig. S2A; see also LaRue et al., 2015; Schoofs et al., 2010). Fore tibia ratios (FT) follow similar (although not the same) relationships relative to D. melanogaster (see green bars Fig. 2B and Fig. S2A), suggesting that fore tibias (for simplicity referred as tibias thereafter) can be used as a proxy for body size. However, the wing ratio (W) in D. virilis remains almost unchanged relative to D. melanogaster (W=1.02, females and W=1.1, males; magenta bar in Fig. 2B and Fig. S2A). We then proceeded to obtain the population scaling relationships (static allometry) between wings and tibias of each species using the SMA model II regression (see Materials and Methods; Table 1, Fig. 2C; Fig. S2B).

To make a comparison across species, we computed the 95% confidence intervals of the slope and intercept of each scaling relationship (Fig. 2D,E). Notably, the slope of D. melanogaster exhibits an isometric behavior (slope of 1 contained in the 95% confidence interval), but D. simulans, D. ananassae and D. virilis all exhibit slopes larger than 1 (Table 1). Although this slope effect is interesting on its own (see Discussion), we decided to focus our attention on a more dramatic result, namely, the difference in intercept between D. virilis and the rest of the species (Fig. 2C,E). Because the intercept is arbitrarily computed at 0, which is far away from the average tibia size, we computed the intercept of the scaling relationships at the median value of the tibia for all the data (dashed line in Fig. 2C and Fig. S2B). Remarkably, the 95% confidence intervals of the intercepts at the median show a dramatic difference in relative wing size.

We considered the possibility that the difference in the intercept of these wing-to-tibia scaling relationships was controlled by cis-regulatory elements of the vg gene (Fig. 1C). From the two intronic cis-regulatory elements that drive vg expression (Fig. 3), we focused on the vg^QE because it controls expression in most of the wing pouch (compare cyan and magenta patterns in Fig. 3A), and previous work has implicated this element in wing growth (Muñoz-Nava et al., 2020; Parker and Struhl, 2020; Zecca and Struhl, 2007b). We first examined the sequence of the reported vg^QE (Dworkin et al., 2009; GenBank ID: FJ513071.1), in order to determine whether a cis-regulatory element also exists in the other three Drosophila species. We then aligned the sequences of the whole fourth intron of each four species selected. The alignment recovers sequence conservation and the presence of regulatory domains to promote vg transcription (as the Scalloped Interaction Domain, SID), an aspect that allowed us to predict a vg^QE in each of these species (Fig. 3B).

To evaluate the influence of the vg^QE in Drosophila wing allometry, we use the CRISPR/Cas9 Homologous Directed Repair (HDR) system (Evans, 2017; Port et al., 2014) to replace the vg^QE sequence in D. melanogaster with the predicted vg^QE from other species. We decided to replace only with the vg^QE of D. virilis because it was the only species in which we obtained a dramatic change in wing-to-tibia allometries (Fig. 2C and Fig. S2B). Thus, we designed a guide RNA (gRNA) expression vector using a pCFD4 backbone (Howard et al., 2021; Port et al., 2014) that contains two gRNA sequences to target the vg^QE region in D. melanogaster (Fig. 4A). We then generated a vg^QE homologous donor plasmid containing the sequence of the predicted vg^QE of D. virilis and D. melanogaster (as a control) to act as a template for the HDR repair system (Fig. 4A). The experimental design of this donor construct allowed us to generate a base vector, which contains left and right D. melanogaster homologous arms and DsRed, that could work as a single vector to produce different vg^QE replacements using the same gRNA sequences and donor backbone of pHD-DsRed (Fig. 4A). Plasmids with the gRNAs and the D. virilis or D. melanogaster vg^QE vectors were verified through PCR (Fig. S3) and injected into embryos that maternally carry the Cas9 endonuclease (nos-Cas9 flies) to generate the vg^QE-replaced flies.

After obtaining stable stocks in which the vg^QE of D. melanogaster (vg^QEDmel) or D. virilis (vg^QEDvir) was replaced with the endogenous vg^QE of D. melanogaster, we asked whether these replacements could affect the wing-to-tibia scaling relationships. With this aim, we generated homozygous stocks (except in the case of vg^QEDmel, for which homozygous animals were not viable and we used vg^QEDmel/+ flies instead) and plotted length traits (Fig. 4B; Fig. S4A) and static scaling relationships as we did for wild-type animals, using the model II regression (SMA) (Fig. 4C–E and Fig. S4B–D). Tibia measurements are only slightly (although significantly) different in vg^QEDmel and vg^QEDvir flies (Fig. 4B; Fig. S4A), but vg^QEDvir animals have a marked relative reduction in wing length with respect to vg^QEDmel controls and wild-type D. melanogaster itself (Fig. 4B; Fig. S4A). This suggests that the vg^QE of D. virilis likely carries genetic information to reduce relative wing size. However, this effect was not as dramatic as in wild-type D. virilis, as revealed by the wing-to-tibia scaling relationships (Fig. 4C; Fig. S4B), suggesting that other inputs are needed to achieve the degree of relative wing size reduction that D. virilis exhibits.

Following the analysis that we performed in wild-type flies, we computed the 95% confidence intervals of slopes and intercepts in the CRISPR-replaced stocks (Fig. 4D,E). Our analysis shows that vg^QEDvir animals have a marked reduction in wing length with respect to vg^QEDmel controls, although they do not exhibit the same effect in the intercept as in wild-type D. virilis flies. Nonetheless, this experiment reveals that a single cis-regulatory element has the potential to modulate scaling relationships in animals.

The evolution of scaling relationships between body parts is a major source of phenotypic diversity in animals. The diversity of wing to whole-body relationships among insect species displays an astonishing repertoire of allometries; butterflies have enormous wings compared with their bodies, while bees have quite the opposite wing-to-body ratio. While the evolution of these scaling relationships is likely to be very complex among distant species, the evolution of wing to whole-body scaling in more closely related species may be pinpointed to the evolution of specific regulatory changes in the genome. The cis-regulatory hypothesis, in which phenotypes evolve through changes in cis-regulatory regions of key genes has received a lot of attention in the literature (Carroll et al., 2005; Jiggins et al., 2017; Rebeiz and Tsiantis, 2017; Wray, 2003, 2007); however, very few examples in which evolutionary changes have been mapped to specific cis-regulatory interactions have been reported (Stern and Orgogozo, 2009; Wittkopp and Kalay, 2011).

In this work, we used CRISPR-mediated replacement (Fig. 4A) of a specific regulatory sequence within the vg wing selector gene from D. virilis, which diverged from D. melanogaster at least 40 million years ago (Powell and DeSalle, 1995) and displays a reduction of wings relative to the rest of the body (Figs 2 and 4), into the genome of D. melanogaster. We found that this cis-regulatory element itself does contribute, at least partially, to a change in the intercept of the allometric relationship (Fig. 4C; Fig. S4B), providing support for the hypothesis that evolution of scaling relationships could be driven by genetic variations in the vg^QE sequence. It is unclear, however, whether the common ancestor of these species had a relative wing size such as in species of the melanogaster group or a smaller relative wing size such as in D. virilis, but our data support that an important evolutionary event at the vg^QE took place after the split of these groups (Fig. 5).

Although we do not focus on changes in slopes in this study, we reported statistically significant differences in the slopes of allometric relationships (Table 1 and Fig. 2C). The elucidation of slopes in wing-to-tibia scaling relationships of these and other Drosophila species might shed light on how traits generate the characteristic diversity of individual species (Wilcox et al., 2023). However, future studies examining slope changes of specific species in response to environmental effects such as temperature or nutrition may reveal how differences in relative wing sizes among species evolved (Shingleton et al., 2009).

How could the replacement of a cis-regulatory element affect the ability to scale the size of an organ with respect to the size of another reference within the animal? Under the experimental conditions evaluated, the vg^QEDvir may adjust its response to Dpp and Wg, morphogens that have been shown to directly affect wing development in D. melanogaster (Parker and Struhl, 2020). The reduced growth of wings in vg^QEDvir mutants could be a result of changes in morphogen spread or degradation in the common ancestor of D. virilis and D. melanogaster. In this context, mathematical modeling approaches underlying the dynamics of morphogen gradients that govern wing development could serve as hypothesis-generation tools to determine the possible evolutionary trajectories driving the observed wing diversity in these Drosophila species.

Another approach for inferring possible molecular mechanisms that could explain the observed changes in relative wing size is through bioinformatic analysis. For instance, in addition to the participation of Wg and Dpp in QE activation and vg autoregulation (Parker and Struhl, 2020; Zecca and Struhl, 2007a), a previous study in D. melanogaster and D. virilis suggest that Dpp has a positive role and Wg has a negative role in the regulation of an enhancer that promotes visceral mesoderm induction (Lee and Frasch, 2005). A DFR (Drosophila Drifter) POU Domain transcription factor has also been described in D. virilis and D. melanogaster that binds to an adjacent site in QE MAD-binding site that may be essential to QE activation in the wing pouch, through a complex between MAD and DFR (Certel et al., 2000). This MAD site is apparently also a binding site for Brinker (Brk), presumably acting as a competing repressor of vg transcription (Kirkpatrick et al., 2001). Perhaps in the vg^QEDvir flies, Brk-mediated repression is somehow favored, altering morphogen spreading and the relative size of the wing. This vg^QEDvir repressor hypothesis could also explain why when we delete the vg^QE sequence in D. melanogaster, the wing phenotype is fully rescued (Farfán-Pira et al., 2022), presumably through the action of shadow enhancers (Hong et al., 2008), but in the vg^QEDvir replacement the wing phenotype is not.

The vg^QE sequence has also been studied to determine which components are necessary to induce vg transcription, where SID domains of Scalloped (Sd) interact with DNA consensus sequence (TGGAATGT) located in the QE region (Halder et al., 1998; Klein and Arias, 1999; Simmonds et al., 1998). As we noted in our alignment (Fig. 3B and Fig. S5), SID domains are evolutionarily conserved in the four species examined, suggesting that this role of Sd to form a complex with Vg is a conserved mechanism to induce the transcriptional mechanism of vg expression in Drosophila. However, additional sites within the QE also influence vg transcription; for example, the Drosophila transcription factor Mothers Against Dpp (MAD) has an N-terminal homology region (mad1) that binds within the QE (GCCGnCGC) sequence in D. melanogaster and has a role in direct expression of vg in the pouch and wing formation. Furthermore, mutation of this consensus sequence and the inhibition of interaction between MAD and the QE sequence produces animals with smaller wings (Certel et al., 2000; Kim et al., 1997). According to a previous report, in D. virilis, the MAD binding site is present in the sequence of regulatory sequences of the daughters against dpp (dad) gene (Weiss et al., 2010), suggesting that D. virilis may have similar mechanisms of vg regulation by MAD in QE sequences, but further examination is necessary to confirm or refute this hypothesis.

Strikingly, the presence of vg^QE sequences located in the fourth intron of the vg gene in the honey bee Apis mellifera (which diverged from Drosophila about 300 million years ago) shows similar patterns of expression in imaginal wing discs (Prasad et al., 2016). Thus, the evolutionary conservation of QE in other insect groups suggests that the cis-regulatory role of the vg^QE in driving wing scaling relationships may be relevant in more distant species as well.

We thank José Luis Fernández-López and Rafael Rodríguez-Muñoz for technical assistance. We are grateful to the core facilities at Cinvestav-IPN and the University of Arkansas for providing us access to equipment to conduct the experiments for this study. We thank Eduardo Ochoa for taking and sharing the photos of the specimens in Fig. 2. We also thank Ehab Abouheif, Fanis Missirlis, Augusto Cesar Poot-Hernández and members of the Nahmad and Evans laboratories for discussions.