The D. melanogaster query genes were retrieved from the compilations of insect developmental genes published by the Tribolium Genome Sequencing Consortium (Supplementary Tables 11, 13 in Richards et al., 2008). D. melanogaster amino acid sequences were retrieved from GenBank. The genome databases used in this study included Drosophila melanogaster genome database version 5.2 (Adams et al., 2000), Anopheles gambiae str. PEST genome database version 2.2 (Sharakhova et al., 2007), Tribolium castaneum Georgia GA2 genome database version 3.0 (Richards et al., 2008), Apis mellifera DH4 genome database version 4.0 (Honeybee-Genome-Sequencing-Consortium, 2006), and Acyrthosiphon pisum genome assembly 1.0 (International Aphid Genomics Consortium, 2010).

The expanded searches for conserved D. melanogaster DGDs in other dipteran genomes were conducted in Mayetiola destructor genome assembly 1.0 (Zhao et al., 2015), Lutzomyia longipalpis genome assembly 0.1 (Sand-Fly-Sequencing-Consortium, 2011), Drosophila virilis genome assembly dvir_caf1 (Drosophila 12 Genomes Consortium et al., 2007), Musca domestica genome assembly MdomA1 (Scott et al., 2014), Stomoxys calcitrans genome assembly ScalU1, Glossina morsitans genome assembly GmorY1 (International Glossina Genome Initiative, 2014), Ceratitis capitata genome assembly Ccap_1.1 (Papanicolaou et al., 2016), Lutzomyia longipalpis assembly LlonJ1 (Sand-Fly-Sequencing-Consortium, 2011), Phlebotomus papatasi genome assembly PpapI1, Aedes aegypti genome assembly AaegL3, and the Rhodnius prolixus genome assembly (RproC3). In the cases where the developmental gene duplication occurred within the Drosophila genus, we searched further species from the Drosophila and Sophophora subgenera (Drosophila 12 Genomes Consortium et al., 2007). All of these ortholog searches for were conducted either in the VectorBase or the NCBI genome databases (Pruitt et al., 2005; Lawson et al., 2009). Complementary BLAST searches were carried out in the Episyrphus balteatus transcriptomes SRX042197, SRX042231, and SRX1131533 (Lemke et al., 2011).

The gene families investigated in this study were defined as monophyletic groups of closely related paralogs in the Drosophila genome, as inferred by a 6-step procedure: 1. Each developmental gene compiled in Richards et al. (2008) served as query seed to collect candidate gene family members by BLASTP (Altschul et al., 1997) against the Drosophila protein sequence database. 2. A maximum e-value of 1.0e⁻¹¹ and a minimal sequence identity D value of 30% were implemented as combinatorial cut-off filter in a first collection of candidate gene family members. 3. Core paralog clusters were extracted from the expansive e-value structured list of candidate gene family members by removing all paralogs below the highest ranked candidate paralog whose the e-value value was smaller than five orders of magnitude than that of the next ranked paralog. 4. To reduce the chance of excluding highly diverged gene family members, the core paralog clusters were retroactively expanded by re-adding the best ranked candidate paralogs from the preliminarily excluded genes until the e-value differed less than five orders of magnitude from the next best hit, indicating saturation of sequence divergence. 5. The gene family membership of each candidate paralog was then assessed by reciprocal BLAST against the Drosophila protein sequence database. Candidate paralogs which returned the Drosophila query seed sequence as top hit were accepted as confirmed gene family members. 6. Candidate gene families with shared members were merged to form non-redundant gene families. This procedure resulted in a total of 377 gene families comprising 661 individual D. melanogaster genes (Supplementary Data File 1).

All members of the Drosophila developmental gene families were used as queries to search the genome databases of mosquito, flour beetle, honeybee, and pea aphid with BLASTP or TBLASTN (Altschul et al., 1997). Putative homologs with an e-value equal or lower than 1.0e⁻⁰⁴ were tested for orthology by reciprocal BLAST against the D. melanogaster RefSeq protein database. Orthology relationships between recovered homologs for a given gene family were further assessed by gene tree analysis. To this end, multiple sequence alignments were generated with ClustalW2 (Larkin et al., 2007) or MUSCLE (Edgar, 2004). Ambiguously aligned positions and divergent regions were removed with Gblocks (Castresana, 2000) at default settings. Tree-Puzzle was used for maximum likelihood tree search (Strimmer and Von Haeseler, 1999; Néron et al., 2009), applying the JTT model of protein sequence evolution and accommodating for rate heterogeneity between sites with four gamma rate categories (Whelan and Goldman, 2001). The majority of these analyses were performed in the now retired Mobyle Project environment (Néron et al., 2009). For a s selection of gene families, maximum likelihood gene trees were generated with MEGA7 (Kumar et al., 2008).

Orthologs of D. melanogaster lineage-specific DGDs in other dipteran species were searched by reciprocal BLAST followed by gene tree analyses. Supplementary Data File 2 contains the sequences of all compiled homologs of Drosophila DGDs.

Non-synonymous (dN) and synonymous substitution (dS) divergences were estimated with the yn00 algorithm of PAML version 3.15 (Yang, 1997). In the case of multiple duplications per gene family, dS and dN of duplicated descendants were averaged.

Relative rate tests were conducted with PHYLTEST 2.0 (Kumar, 1996), applying the Benjamini & Hochberg False Discovery Rate (FDR) correction (Benjamini and Hochberg, 1995) and using singleton homologs from T. castaneum or A. mellifera for outgroup comparison.

Information on gene function was retrieved from FlyBase and the primary literature (Tweedie et al., 2009). Expression patterns were explored in literature compiled through FlyBase and by examining gene specific entries in the FlyExpress image database when available (Kumar et al., 2011).

To explore the impact of gene duplication on the genetic architecture of Drosophila development compared to those of other insects, we explored the duplication histories of 377 conserved developmental gene families. These were represented by 661, 642, 622, 620, and 696 individual genes in Drosophila, Anopheles, Tribolium, Apis, and Acyrthosiphon, respectively (Supplementary Data Files 1, 2). For ~10% of the investigated gene families, reciprocal BLAST results produced evidence of duplications in more than one lineage. In these cases, the phylogenetic relationships between homologs were further examined by gene tree estimation and analysis.

Consistent with the previously reported overall genome duplicate richness of the pea aphid (International Aphid Genomics Consortium, 2010; Shigenobu et al., 2010), the highest number of lineage-specific DGDs was found in the pea aphid, where a total of 93 gene duplications were distributed over 61 gene families (Figure 1). More surprisingly, Drosophila stood out with the second highest number of lineage-specific duplicates, estimated at 62 duplication events in 49 gene families. Considerably fewer lineage-specific DGDs were detected in the remaining three species with 20 duplications in 14 gene families in the honeybee, 22 duplications in 20 gene families in the red flour beetle, and 26 duplications in 21 gene families of the mosquito (Figure 1, Supplementary Data Files 1, 3). Taken together, these results revealed an exceptionally high number of DGDs not only in the pea aphid, but also in the lineage to D. melanogaster.

Previous studies have shown that tandem gene duplication is the major contributor of duplicated genes (~80% of evolutionarily very young duplicates) in Drosophila species followed by retrotransposition (~10%) (Zhou et al., 2008). Consistent with this, all of the Drosophila lineage-specific gene duplicates identified in our previous study of vision-related genes represented tandem duplicated paralogs (Bao and Friedrich, 2009). To further probe the generality of these findings, we explored the frequency of physical linkage among the DGDs sampled from Drosophila, Anopheles, Tribolium, and the honeybee. The pea aphid was not included in this analysis due to the preliminary state of chromosome scaffolds at the time of analysis. In the examined species, over 65% of the sampled sister paralogs were on the same contig. Moreover, between 40 and 60% of DGDs, depending on the species, were physically linked within less than 500 kb (Figure 2, Supplementary Data File 4). Factoring in the expected breakdown of physical linkage over time, these numbers identified tandem gene duplication as the generally predominant source of gene duplicates in insects.

To gain insight into the time course of DGD accumulation in the five examined insect lineages, we calculated evolutionary distances at synonymous sites (dS) between sister duplicates as proxies of DGD ages (Lynch and Conery, 2000) (Supplementary Data File 5). dS distributions were compared after binning into 7 age classes (Figure 3). Again consistent with previous studies (International Aphid Genomics Consortium, 2010), there was a marked peak of gene duplicates in the youngest age class (0 < dS < 1) for the pea aphid, amounting to 93 duplications (55%). This number was seven times higher than the maximum number of DGDs in this age class in any of the other species (A. mellifera: 7) and at least 2.5 times higher than the maximum number of duplications in any other age class across species. Of note, the second highest number in the youngest dS duplicate age class was detected in the honeybee. This, however, was largely due to six rounds of gene duplication in a single gene family (farnesyl pyrophosphate synthases) (Supplementary Data File 1), thus not reflecting a broader trend. Further, consistent with the predicted outcomes of birth-death models of gene duplicate evolution (Lynch and Conery, 2000), the pea aphid gene duplicate number dropped to 15 in the next oldest age class (1 ≤ dS < 2), followed by a milder but consistent decrease over the remaining older age classes, except for a mild secondary peak in the 4 ≤ dS < 5 age bin.

Contrasting with the pea aphid DGD age profile, DGD ages did not peak in the youngest gene duplicate group for any of the other four species. Instead, the numbers of Drosophila, Tribolium, and Anopheles DGDs in older age classes invariably exceeded that in the 0 < dS < 1 class (Figure 3). This trend was most pronounced in Drosophila where the majority of DGDs (88%) were captured in the age class range 2 ≤ dS < 5. Moreover, the number of Drosophila DGDs in this age range exceeded that of any other species, including the pea aphid. This finding suggested that the exceptional number of the Drosophila lineage-specific DGDs had accumulated substantially deeper back in time than in the pea aphid and, as a corollary, represented distinctly more long-term preserved DGDs.

The accuracy of dS divergences as proxies of gene duplicate age decreases with time depth due to substitution saturation and limited sequence sample size in terms of alignable conserved sequence regions. Therefore, to scrutinize the antiquity of the Drosophila DGDs further, we investigated their conservation in nine additional dipteran genomes by reciprocal BLAST searches and gene tree analysis (Supplementary Data File 6). This approach sorted the D. melanogaster DGDs into five age groups (Figure 3):

To explore the functional impact of DGD accumulation in the Drosophila lineage, we capitalized on the rich documentation of gene function in Drosophila, which allowed us to parse the Drosophila lineage DGDs into three functionalization groups (Supplementary Data File 8): (I) Fully redundant defined by the lack of detectable differences in spatiotemporal expression between sister paralogs, restriction of detectable phenotypic abnormalities to animals doubly mutant for both gene family members, or both; (II) Partially redundant defined by partial overlap of expression patterns, phenotypic abnormalities unique to Markus both double-mutant and paralog-specific mutant animals, or both; (III) Functionally independent as evidenced by the lack of overlapping expression patterns, lack of evidence of compensatory genetic interactions in the literature, or both.

Sufficient information on gene expression, function or both was accessible for 48 of the Drosophila lineage DGD paralog pairs, representing 37 gene families due to multiple duplications in 10 gene families (Supplementary Data File 8). Of these, 4 (8%) were characterized as fully redundant, 25 (50%) as partially redundant, and 20 (42%) as fully non-redundant. Proportions, however, varied when the categories were parsed by gene family age groups (Figure 5).

All three functionally characterized maximally 30 million years old DGDs were non-redundant. In the 10 functionally characterized 30–65 million years old DGDs, however, only 6 were documented as non-redundant while 4 were documented as partially redundant and 1 as fully redundant. The proportion of fully non-redundant paralogs was even lower in the 65–80 million years age group with only 2 non-redundant paralogs, compared to 8 redundant paralogs. The large 80–230 million years age group contained 8 non-redundant paralogs, 8 partially redundant paralogs, 2 partially redundant paralog triplets, and, most notably, 3 fully non-redundant paralog pairs. Finally, one example each of partial redundancy and non-redundancy was found in the group of 230–240 million years old DGDs (Figure 5).

Overall, thus, functional evidence from Drosophila genetics suggests a substantial amount of long-term conserved genetic redundancy in the Drosophila DGDs, which did not decline over time. Instead, a substantial and consistent number of the Drosophila lineage DGDs maintained their likely ancestral genetic redundancy for up to 200 million years, resulting in an approximate balance of diverged vs. redundant DGD paralog fates.

Finally, to gauge the impact of non-redundant DGDs in the Drosophila lineage, we determined the proportion of significantly asymmetrically sequence diverged sister paralog pairs in the Drosophila lineage DGDs via relative rate tests (Supplementary Data File 9). Following gene duplication, loss of genetic redundancy may occur due to complete subfunctionalization with little or no phenotypic consequences, and hence a conservative, trajectory, or neofunctionalization, a gene regulatory and potentially phenotypically innovative trajectory. As a rule of thumb, neofunctionalization is often associated with a transient, yet dramatic and significant, acceleration of protein sequence change compared to the ancestrally functioning paralog. Relative rate tests therefore serve as an efficient approach to identify candidate neofunctionalized paralogs (Conant and Wagner, 2003).

After correcting for multiple testing, close to 50% of all Drosophila DGD paralog pairs were found significantly asymmetrically diverged. The same was true for the subcohort of functionally characterized DGDs, in which case 24 out of 49 were diagnosed to have significantly asymmetrically diverged. As expected, non-redundant duplicates were characterized by the highest proportion of asymmetrically diversified sister paralogs, i.e., 60%, followed by partially redundant paralogs with 44%. Only one of the four fully redundant sister paralogs was marginally significantly asymmetrically diverged (Supplementary Data File 9).

Parsed by age groups, the proportion of significantly asymmetrically diverged DGDs varied from 30% (80–230 million years old) to 100% (230–240 million years old) between DGD age groups (Figure 5). The most strongly diverged DGD paralogs, however, were contained in the youngest age groups of 0–35 and 30–65 million years old DGDs (Figure 5). The analysis of relative sequence divergence between Drosophila DGD sister paralogs thus uncovered tentative evidence of a higher rate of neofunctionalization, and hence phenotypically innovative DGD trajectories, in the past ~60 million years of Drosophila lineage evolution, contrasting with the pronounced degree of functional redundancy among the more ancient Drosophila lineage DGDs reaching back to up to 200 million years.

For the over 350 developmental gene families investigated in this study, Drosophila and the pea aphid stand out with substantially higher percentages of lineage-specific duplications (13.5 and 16.2%, respectively) compared to Anopheles (5.57%), Tribolium (5.31%), and Apis (3.71%). As noted, this result is consistent with the known genome-wide preponderance of duplicated genes in the pea aphid (International Aphid Genomics Consortium, 2010). For Drosophila, however, previous genome-wide studies did not report evidence of notable differences in duplicate numbers compared to other insect genome models (Zdobnov et al., 2002; Honeybee-Genome-Sequencing-Consortium, 2006; Richards et al., 2008; International Aphid Genomics Consortium, 2010). One possible explanation is that this difference is specific for developmental genes and not a genome-wide phenomenon in Drosophila and related Diptera. A notably higher number of duplicated genes, however, has also been found for structural vision genes (Bao and Friedrich, 2009), raising the possibility of a more general scope. Genome-wide surveys of lineage-specific gene duplication will provide an ultimate answer to this question.

Further confidence in the accuracy of our comparative gene family analysis comes from consistent findings in earlier, gene-specific studies. In total, 14 (~30%) of the Drosophila lineage-specific expanded developmental gene families covered here have been previously identified as such (Supplementary Data File 6). Also the taxonomic distribution of recently identified whole genome duplication events in the insect tree of life is consistent with our findings (Li et al., 2018).

While the overall evidence is compelling that DGDs played an exceptional role during the evolution of brachyceran Diptera, resolving its timeline to a satisfactory degree will require substantial further work. Based on the relatively high number of DGDs associated with the basal-most branch in the schizophoran Diptera covered in our analysis (Figure 4), it is, for instance, tempting to speculate that DGD accumulation may have spiked in conjunction with the origin of cyclorrhaphan Diptera (Figure 6). This inference, however, will require analyses of DGD conservation in the cyclorrhaphan family cluster Platypezoidea and ancient cyclorrhaphan key families such as the Syrphidae (hoverflies) (Figure 6).

In preliminary studies, we searched embryonic and adult transcriptome data of the hoverfly Episyrphus balteatus (Supplementary Data File 5) (Lemke et al., 2011). The results from this exercise suggest that at least 8 of the 22 gene families in the 80–230 million years time window duplicated prior to the diversification of the Cyclorrhapha, implying that 14 families might have expanded specifically in the cyclorrhaphan stem lineage after its separation from the ancestor of modern Syrphidae (Figure 6). However, in the absence of whole genome coverage, the latter number may be an overestimate and we therefore abstained from including these findings in our current gene duplication tree (Figure 4).

The exceptional accumulation of DGDs in the brachyceran clade of the Dipteran tree of life prompts questions regarding its impact on the genetic control of development and, by extension, body plan evolution. The proximate effect to be expected from DGD accumulation is the emergence of novel gene regulatory network components, a prediction that has been documented for select Drosophila lineage DGDs. The cyclorrhaphan-specific Hox3 transcription factor paralog bicoid (bcd), for instance, is a paradigm example of an extremely asymmetrically evolved, neofunctionalized DGD. As novel regulator of early anterior patterning in the Drosophila embryo, bcd interacts with a rich array of ancient, pre-dipteran regulators. These interactions include the RNA-binding protein Exuperantia, which predates Brachycera and insects (MacDonald et al., 1995; de Oliveira et al., 2017) and direct target genes as ancient as Orthodenticle (Finklstein and Perrimon, 1990), hunchback (Driever and Nüsslein-Volhard, 1989; Finklstein and Perrimon, 1990), and caudal (Wolff et al., 1998).

The emergence of new gene regulatory network components in turn is predicted to facilitate, or to be driven by, advantageous phenotypic change. In the case of developmental regulators, this can come in the form of changes in body plan traits or in their development. Likewise consistent with this prediction, DGD-associated patterning innovations are well-documented for brachyceran DGDs. The neofunctionalization of bcd, as a case in point, occurred in the context of the dramatic compaction of two ancestral extraembryonic membranes, amnion and serosa, into a single one, the amnioserosa, in the lineage to cyclorrhaphan Diptera (Stauber et al., 1999; Rafiqi et al., 2008). Similarly timed expansions of signaling pathway-related gene families likewise contributed to the regulatory evolution of the amnioserosa (Richards et al., 2008; Fritsch et al., 2010; Lemke et al., 2011). The expansion of the achaete-scute complex, which predates the schizophoran radiation (Negre and Simpson, 2009), affected the evolution of thoracic bristle patterns (Skaer et al., 2002). The same holds for the expansion of the Drosophila lineage-specific expansion of the enhancer of split gene complex (not included in this analysis) (Baker et al., 2011).

Altogether, our analyses identified close to 30 brachyceran DGDs with asymmetrically diverged protein sequences, which thus potentially produced novel functionalities. While even dramatically asymmetrically diverged DGDs can maintain genetically redundant ancestral functions (Bao et al., 2012), it is reasonable to conclude from the large number of asymmetrically diverged DGDs that developmental gene family expansions did play a innovative roles at specific stages of brachyceran body plan evolution. Our data further indicate a higher proportion of phenotypically innovative DGD trajectories in the past ~60 million years of the brachyceran lineage leading to Drosophila (Figure 4). This combines with tentative evidence of peaking DGD accumulation in the basal Cyclorrhapha and Schizorrhapha, which gave rise to over 40% of extant dipteran diversity (Yeates and Wiegmann, 1999, 2007) and acquired an exceptionally large number of body plan changes (McAlpine, 1989; Lambkin et al., 2013) (Figure 6). The same can be stated for two younger expansive subclades nested within the Cyclorrhapha: the Schizophora and the Calyptratae (Figure 6). The exceptional accumulation of DGDs preceding the diversification of cyclorrhaphan Diptera may thus be functionally related to the dramatic changes of this clade with regards to embryonic and postembryonic development as well as overall developmental speed. Intriguingly, an acceleration of development could at the same time explain the long-term conserved genetic redundancy of many Drosophila DGDs.

Besides likely or known neo- or subfunctionalization DGD trajectories, our analyses uncovered a fairly balanced mix of partially or fully redundant vs. functionally completely diverged Drosophila DGDs despite their overall antiquity. As a case in point, one of the two oldest Drosophila DGDs sampled in this study, the zinc finger gene sister paralog pair disco and discor, constitutes an exceptionally well-studied example of genetic redundancy (Heilig et al., 1991; Mahaffey et al., 2001). The Drosophila disco paralog has been found to function in the developing visual system (Steller et al., 1987; Lee et al., 1991; Campos et al., 1995), early embryonic segment identity specification (Robertson et al., 2004), and leg development (Dey et al., 2009). Comparative analyses of the disco/discor singleton ortholog in Tribolium have led to the conclusion that the leg and visual system patterning functions of Drosophila disco are ancestral for higher insects while the embryonic segment identity specification originated at a later point in time (Patel et al., 2007). In the context of embryonic gnathal head segment development, disco and discor are fully redundant (Mahaffey et al., 2001). In the context of leg development, disco and discor are hypothesized to be partially redundant (Dey et al., 2009), consistent with their precisely overlapping expression patterns in the leg imaginal disks (Mahaffey et al., 2001). Of further note is the apparently conserved linkage of the two genes, which are separated by less than 100 kb on the Drosophila X-chromosome (Mahaffey et al., 2001) and less than 200 kb on scaffold NW_004523853 in C. capitata. Combined with the conservation of both paralogs in the Hessian fly (Figure 4), these data point at potentially over 200 million years of preserved redundant regulation of head segmentation and leg patterning by disco and discor in brachyceran Diptera.

The large number of early Brachycera-specific DGDs includes three additional examples of long-term conserved genetic redundancy: The homeobox transcription factor duo Bar-H1 and Bar-H2 (Higashijima et al., 1992a,b), the zinc finger transcription factor pair spalt major and spalt related (Barrio et al., 1999; Elstob et al., 2001; Cantera et al., 2002; Dong et al., 2003), and the Dorsocross 1-3 T-box transcription factor paralogs (Reim et al., 2003). These fully redundant DGD paralogs are joined by 12 partially redundant paralogs that originated prior to the diversification of the Schizophoran clade (Figures 4, 5). Moreover, partially redundant paralogs continue to represent the majority of DGDs that originated during early schizophoran diversification 65–80 million years ago. They also represent a considerable fraction of the DGDs that originated before the origin of drosophilid Diptera (Figure 4). Even though the exact proportion of genetically redundant vs. non-redundant interactions is likely overestimated due to the usually focused and therefore inherently incomplete nature of gene function studies, the substantial proportion of long-term conserved genetic redundancy in the brachyceran DGDs raises the question of possible adaptive aspects.

The fitness benefits of long-term conserved genetic redundancy have been studied for considerable time (Krakauer and Nowak, 1999; Bessa et al., 2009; Payne and Wagner, 2015). Neutral models predict that gene paralogs eventually diverge by differential loss of functionalities, which has found support in large-scale analyses of expression domain evolution in duplicated genes (Lynch and Conery, 2000; Oakley et al., 2006; Mendonca et al., 2011). Recent gene-specific and genome-wide studies, however, produced evidence for a role of genetic redundancy in securing developmental, genetic, and environmental robustness over hundreds of millions of years in part through conservation of genetically redundant gene paralogs (Celniker et al., 2002; Maslov et al., 2004; Pasek et al., 2006; Dean et al., 2008; Vavouri et al., 2008; Yang et al., 2009; Bao et al., 2012; Buscà et al., 2015). As an example, the recent discovery of partially, yet long-term conserved redundant roles of paralogs of the Drosophila lineage expanded MADF-BESS transcription factor family in the development of the wing hinge has been proposed to be explained by the benefit of developmental robustness (Shukla et al., 2014).

An attractive explanation for the apparent gene duplication facilitated increase of genetic robustness in brachyceran lineages is the previously noted trend toward increased developmental speed in the higher Diptera, as reflected, for instance, by the transformation from short to long germband development during early cyclorrhaphan evolution (Tautz et al., 1994). An acceleration of complex pattern formation processes can be envisioned to impose increased demands on regulatory precision and robustness. Further consistent with the notion of accelerated development in the higher Diptera are the aforementioned compaction of extraembryonic membranes and the prevalent expression of short, intron-free transcripts during early embryonic development in Drosophila (Nunes da Fonseca et al., 2010). As a specific example, the zinc finger transcription factor paralogs knirps (kni) and knirps-related gene (knrl) are coexpressed in the early Drosophila embryo. Remarkably, however, only kni can execute the correlated gap gene patterning function because the ~20-fold longer intronic sequence of knrl prevents the on-time completion of transcript formation imposed by the fast cell cycle succession during early Drosophila embryogenesis (Rothe et al., 1992; Swinburne and Silver, 2008). In this light, it becomes tempting to speculate that the widespread presence of redundant enhancer elements may represent a lineage-specific corollary of accelerated development in the higher Diptera (Hong et al., 2008; Frankel et al., 2010; Perry et al., 2010; Wunderlich et al., 2016). From a practical point of view, this conjecture predicts that other genomic models of insect development might less replete with redundant gene functionality compared to Drosophila, which would be good news for ongoing large scale gene knockdown projects in non-dipteran insect species such as Tribolium (Dönitz et al., 2015).

Finally, while providing developmental robustness as a proximate benefit, genetic redundancy has been found to serve as critical source for new gene regulatory opportunities in species diversification and body plan evolution (Wagner, 2008; Melzer and Theißen, 2016; Wei and Zhang, 2017). It therefore seems reasonable to hypothesize that the interplay of gene duplication, developmental robustness, and adaptive opportunities played an important role during the vast diversification of brachyceran Diptera.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.