ORIGINAL RESEARCH published: 09 October 2020 doi: 10.3389/fgene.2020.565836 The Cyclically Seasonal Drosophila subobscura Inversion O7 Originated From Fragile Genomic Sites and Relocated Immunity and Metabolic Genes Charikleia Karageorgiou* , Rosa Tarrío* and Francisco Rodríguez-Trelles* Grup de Genòmica, Bioinformàtica i Biologia Evolutiva (GGBE), Departament de Genètica i de Microbiologia, Universitat Autonòma de Barcelona, Barcelona, Spain Chromosome inversions are important contributors to standing genetic variation in Edited by: Drosophila subobscura. Presently, the species is experiencing a rapid replacement of Pedro Simões, high-latitude by low-latitude inversions associated with global warming. Yet not all low- University of Lisbon, Portugal latitude inversions are correlated with the ongoing warming trend. This is particularly Reviewed by: Stephen Wade Schaeffer, unexpected in the case of O7 because it shows a regular seasonal cycle that peaks in Pennsylvania State University (PSU), summer and rose with a heatwave. The inconsistent behavior of O7 across components United States Priscilla Erickson, of the ambient temperature suggests that is causally more complex than simply due to University of Virginia, United States temperature alone. In order to understand the dynamics of O7 , high-quality genomic Emma Berdan, data are needed to determine both the breakpoints and the genetic content. To fill University of Gothenburg, Sweden this gap, here we generated a PacBio long read-based chromosome-scale genome *Correspondence: Charikleia Karageorgiou assembly, from a highly homozygous line made isogenic for an O3+4+7 chromosome. charikleia.karageorgiou@uab.cat Then we isolated the complete continuous sequence of O7 by conserved synteny Rosa Tarrío rosamaria.tarrio@uab.cat analysis with the available reference genome. Main findings include the following: (i) Francisco Rodríguez-Trelles the assembled O7 inversion stretches 9.936 Mb, containing > 1,000 annotated genes; franciscojose.rodrigueztrelles@uab.cat (ii) O7 had a complex origin, involving multiple breaks associated with non-B DNA- Specialty section: forming motifs, formation of a microinversion, and ectopic repair in trans with the This article was submitted to two homologous chromosomes; (iii) the O7 breakpoints carry a pre-inversion record of Evolutionary and Population Genetics, fragility, including a sequence insertion, and transposition with later inverted duplication a section of the journal Frontiers in Genetics of an Attacin immunity gene; and (iv) the O7 inversion relocated the major insulin Received: 26 May 2020 signaling forkhead box subgroup O (foxo) gene in tight linkage with its antagonistic Accepted: 09 September 2020 regulatory partner serine/threonine–protein kinase B (Akt1) and disrupted concerted Published: 09 October 2020 evolution of the two inverted Attacin duplicates, reattaching them to dFOXO metabolic Citation: Karageorgiou C, Tarrío R and enhancers. Our findings suggest that O7 exerts antagonistic pleiotropic effects on Rodríguez-Trelles F (2020) The reproduction and immunity, setting a framework to understand its relationship with Cyclically Seasonal climate change. Furthermore, they are relevant for fragility in genome rearrangement Drosophila subobscura Inversion O7 Originated From Fragile Genomic evolution and for current views on the contribution of breakage versus repair in shaping Sites and Relocated Immunity inversion-breakpoint junctions. and Metabolic Genes. Front. Genet. 11:565836. Keywords: non-B DNA, genome fragility, foxo (forkhead box subgroup O), Akt1 (serine/threonine–protein kinase doi: 10.3389/fgene.2020.565836 B), Attacin antibacterial genes, immunometabolism, thermal adaptation, seasonal selection Frontiers in Genetics | www.frontiersin.org 1 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance INTRODUCTION 2019). The NEHJ mechanism rests upon the occurrence of two or more DSBs. But the source of the DSBs (whether environmental, Chromosome inversions are arguably the genetic traits with the such as ionizing radiation, or spontaneous, such as non-B earliest and richest record of associations with climate (Hoffmann DNA-associated sequence instability, where non-B DNA denotes and Rieseberg, 2008). Research into evolutionary responses to any DNA conformation that is not in the canonical right-handed contemporary global warming (Hughes, 2000; Parmesan, 2006) B form; Lobachev et al., 2007; Zhao et al., 2010; Farré et al., 2015), is therefore faced with the challenge of understanding how the relative contributions of breakage versus repair to shaping inversions originate and spread in populations (Kirkpatrick, breakpoint junctions (Ranz et al., 2007; Kramara et al., 2018; 2010), while trying to determine their roles in climatic adaptation Scully et al., 2019), and the relative frequency with which the (Gienapp et al., 2008; Messer et al., 2016). joined broken ends are from the same chromatid (isochromatid Chromosome inversions are ubiquitous chromosomal model) versus two distinct sisters (chromatid model) (Ranz et al., mutations consisting in the reversal of the orientation of a 2007) or even, as has been more recently suggested by Orengo chromosome segment. They originate through either of two et al. (2019), non-sister chromatids (chromosome model) are major mechanisms, each with its associated distinctive footprints. additional open questions. The first mechanism is intrachromatidal non-allelic homologous Inversions can have direct or/and indirect functional effects recombination (NAHR) between inversely oriented repeats. (Kirkpatrick, 2010). Direct effects are those ascribable to the This mechanism generates inversions with duplications at their mutation per se, as it altered the structure or expression of ends in both the inverted and uninverted states (Cáceres et al., functional sequences at the breakpoints, or the functional 1999). The second mechanism is chromosomal breakage and neighborhood of genes in the cell nucleus (McBroome ectopic repair via non-homologous end joining (NHEJ). This et al., 2020). Indirect effects emanate from their associated mechanism either does not generate duplications or generates recombination–suppression effects when in heterozygous them but at the ends of the inverted state only. These two condition, whereby they can bind together into close linkage types of NHEJ footprints have been explained in terms of association particular combinations of alleles at genetically differences in the mode of breakage. Two modes of breakage distant loci. The evolutionary significance of polymorphic have been advanced: “cut-and-paste” via clean double-strand inversions is often thought to chiefly stem from their indirect breaks (DSBs) that generate blunt ends and staggered. NHEJ effects (Dobzhansky, 1947; Wasserman, 1968; Kirkpatrick and inversions without duplications at their ends would originate Barton, 2006). Although data have been lacking on the relative via cut-and-paste (Wesley and Eanes, 1994), whereas those importance of the two types of effects, there has been renewed with inverted duplications at their ends would originate via interest in using genomics to determine mechanisms for the staggered breaks in one or the two breakpoints. Two staggering spread, establishment, and maintenance or fixation of inversions models for the origin of the inverted duplications have been (Corbett-Detig and Hartl, 2012; Corbett-Detig, 2016; Fuller proposed (Kehrer-Sawatzki et al., 2005; Matzkin et al., 2005; et al., 2016, 2017, 2019; Cheng et al., 2018; Said et al., 2018; Ranz et al., 2007): according to the isochromatid model, the Lowry et al., 2019). Because they usually involve many genes, duplications would be the filled-in single-stranded overhangs chromosome inversions have enhanced potential for affecting that would result from paired single strand breaks (SSBs) located multiple traits, which should expand the opportunities for their staggered with each other on opposite strands of the same maintenance via balancing selection. The extent to which that is chromatid (Kehrer-Sawatzki et al., 2005), whereas according to the case and the types and transience of the balancing selection the chromatid model, the duplications would result from unequal mechanisms involved are only beginning to be elucidated (Kapun exchange between paired sister chromatids, each with one of and Flatt, 2018; Wellenreuther and Bernatchez, 2018; Faria et al., two paired staggered DSBs at each breakpoint (Matzkin et al., 2019). Amid these unknowns, the inversion polymorphisms of 2005). Note that here the terms isochromatid and chromatid have Drosophila subobscura emerged among the first genetic traits switched meanings relative to how they are used in cytogenetics identified as involved in a species’ adaptation to contemporary (Savage, 1976). The two staggering models are chromatid models climate warming (Rodríguez-Trelles and Rodríguez, 1998, 2007; because they assume that inversions originate from either single Balanyà et al., 2006; Rezende et al., 2010). chromatids during premeiotic mitosis (isochromatid), or paired Drosophila subobscura is a native Palearctic species broadly sister chromatids from the same chromosome during meiotic distributed in Europe and the newly invaded areas of North and prophase (chromatid) (Ranz et al., 2007). The models cannot be South America (reviewed in Krimbas, 1992), where it is found distinguished based on the pattern of inverted duplications. Yet generally associated with woodland habitats. It belongs in the the chromatid model has been favored over the isochromatid obscura group, within which it clusters with the recently derived model, because of the length of DNA that would need to be small-island endemics Drosophila guanche and Drosophila unwound by enzymatic activity in the latter model (Ranz et al., madeirensis, forming the subobscura three-species subgroup 2007). The chromatid model is also not without potential (Bächli, 2020). D. subobscura has one of the smallest and least caveats because NHEJ was found to be suppressed during the repetitive Drosophila reference genomes obtained thus far, which meiotic prophase in Drosophila (Joyce et al., 2012; Hughes et al., is distributed among five large telocentric chromosomes (A, J, 2018). The prevalence and distribution of the NAHR and NEHJ U, E, and O) and one small dot (Karageorgiou et al., 2019). mechanisms of inversion formation within and across lineages In stark contrast with its two insular relatives, the species has are currently under debate (Ranz et al., 2007; Delprat et al., evolved highly rearranged chromosome sequences, which is due Frontiers in Genetics | www.frontiersin.org 2 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance to having experienced accelerated fixation rates of paracentric temporal variation; Rodríguez-Trelles et al., 1996) that peaks inversions, especially the A sex chromosome. This situation in summer and drops in winter (Figure 1C). In 2011, it rose has been interpreted as indicative of the inversions’ role in to summer-like levels in spring during a heatwave, with the binding together adaptive alleles in the face of the species’ intense magnitude of the increase closely matching that of the thermal continent-wide gene flow (Karageorgiou et al., 2019). Presently, anomaly (Figure 1C; Rodríguez-Trelles et al., 2013). However, D. subobscura harbors a rich inversion polymorphism, with its (i) the average annual frequency of O7 in northwest Spain five major chromosomes showing parallel adaptive variation remains unchanged after decades of sustained climate warming patterns across latitude (Ayala et al., 1989), seasons (Rodríguez- experienced by the region (Rodríguez-Trelles et al., 2013; our Trelles et al., 1996, 2013), and through a heatwave (Rodríguez- unpublished records). (ii) Following the 2011 heatwave, the Trelles et al., 2013), while rapidly shifting in close association with inversion reached summer-like frequencies in April, but did the ongoing rise in global temperatures (Rodríguez-Trelles and not continue rising through the ensuing summer (Figure 1C), Rodríguez, 1998, 2010; Balanyà et al., 2006). Laboratory attempts perhaps hampered by recessive deleterious alleles (Rodríguez- to establish the causal nature of this association have, however, Trelles et al., 2013). (iii) The Palearctic distribution of O7 largely been inconclusive (Santos et al., 2005; Fragata et al., 2014). is disjointed between the peninsulas of Iberia and Turkey Ultimately, a complete understanding of the role of inversions in (Götz, 1967). These are similar latitude areas separated by adaptation to contemporary climate warming in D. subobscura ∼2,500 km within the continuous species’ range. Assuming will necessarily include the identities and functional properties that the inversion is molecularly the same in the two areas, of the genome sequences affected by them. Advances along this this spatial pattern can hardly be explained on the sole basis line include the isolation and characterization of breakpoint of a postglacial expansion scenario (Menozzi and Krimbas, sequences for 11 of the more than 65 large cytologically visible 1992), considering how rapidly it spread through the recently inversions known for the species, including A2 (Puerma et al., invaded areas of the New World (Prevosti et al., 1988). And 2017), O3 (Papaceit et al., 2012), O4 and O8 (Puerma et al., (iv) in the more studied Iberian Peninsula, the distribution 2016a), E1 and E2 (Puerma et al., 2014), E3 and E9 (Orengo of the inversion has negative or no correlations with the et al., 2015), E12 (Puerma et al., 2016b), and U1 and U2 geographical variation in temperature. For example, the average (Karageorgiou et al., 2019). An overall conclusion is that none annual frequency of the inversion declines from ∼50% to of these inversion breakpoints disrupted any obvious candidate near-zero values along the > 1,000-km stretching from the gene for direct adaptation to temperature, despite the fact northwestern-most to the northeastern-most territories, despite that all but the E3 inversion are supposed to be involved in the latter having a warmer climate than the former (de Frutos, adaptation to climate (e.g., Menozzi and Krimbas, 1992; Rego 1972; Solé et al., 2002; Rodríguez-Trelles et al., 2013). The et al., 2010; Arenas et al., 2018). Apart from the fact that same is true for the West Atlantic fringe of the peninsula thermal traits are genetically complex and that many of the genes along which the inversion levels remain basically the same that impinge on them are still unknown, the above conclusion despite the fact that it stretches seven latitudinal degrees supports that those inversions’ role in thermal adaptation would of steep thermal gradient (Brehm and Krimbas, 1988; Solé be through either position effects, indirect linkage generation et al., 2002; Rodríguez-Trelles et al., 2013). The inconsistent effects, or both. patterns of O7 between components of the ambient temperature As part of a wider effort to develop a high-quality suggest that it is influenced by selective factors other than reference genome for D. subobscura encompassing the species’ temperature alone. rich chromosomal polymorphisms, here we focus on the O7 The O chromosome offers the methodological advantage inversion. The breakpoints of this inversion were located over the other D. subobscura chromosomes that there is an cytologically at subsections 77B/C and 85E on the Kunze– available balancer-strain called Varicose/Bare (Va/Ba) (Sperlich Mühl and Müller standard map (Figure 1A; Kunze-Mühl et al., 1977). In this study, we first used the Va/Ba strain to and Müller, 1958; Götz, 1965). O7 is among the top 10% develop an isogenic line with two identical copies of a wild known largest D. subobscura inversions, stretching most O chromosome carrying the O3+4+7 arrangement. Second, we of the centromere-proximal half of the O chromosome used PacBio long-read technology to generate a high-quality (Krimbas, 1992). In nature, it attains significant frequencies annotated chromosome-scale genome sequence for the line. only in combination with the non-overlapping centromere- Third, we isolated the complete continuous nucleotide sequence distal complex of two overlapping inversions O3+4 , forming of the inversion O7 by conserved synteny analysis of the the chromosome arrangement O3+4+7 (Figure 1B). The tight obtained O3+4+7 chromosome with the available O chromosome association between O7 and O3+4 is likely maintained by from the species’ reference genome, which is structurally O3+4 an interaction between selection and the strongly reduced (Karageorgiou et al., 2019). In addition, we also considered recombination between them (Pegueroles et al., 2010a). two other published sequences of the O chromosome, including O7 could be initially classified as a warm-climate inversion. a high-quality long-read–based sequence from D. subobscura In the Palearctic, it shows a southern distribution. In northwest (Bracewell et al., 2019), and an Illumina-based sequence from Spain, where it has been longitudinally monitored starting in D. guanche (Puerma et al., 2018). We give an account of O7 main mid-1970s (Fontdevila et al., 1983; Rodríguez-Trelles et al., features, together with a detailed description of its mechanism of 1996, 2013), it shows a pronounced regular seasonal cycle formation. Our findings provide clues to the mixed evidence for (estimated to account for more than 60% of the inversion’s this inversion’s role in thermal adaptation. Frontiers in Genetics | www.frontiersin.org 3 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance FIGURE 1 | O7 inversion and O3+4+7 chromosome arrangement of D. subobscura. (A) Light micrograph (400 × ) of the O7 diagnostic loop from two paired polytene O chromosomes of a O3+4+7 /O3+4 heterokaryotype, with indicated cytological map positions of the two inversion breakpoints (Kunze-Mühl and Müller, 1958; Götz, 1965). C and T denote centromere and telomere, respectively. (B) Phylogeny and chromosomal locations of the inversions forming the O3+4+7 arrangement in the subobscura subgroup. Names at the root and tips (bold black) and on branches (bold gray) denote chromosome arrangements and inversions, respectively. The ancestral O arrangement of the subgroup is Oa (Karageorgiou et al., 2019). The chromosome-central inversion Oms (diagonally hatched) is so-called because it became fixed in the last common ancestor of D. madeirensis and D. subobscura (Karageorgiou et al., 2019). In D. subobscura, O3 (blue) and O4 (orange) are two centromere-distal inversions with overlapping cytological map positions originated independently on separate Oms branches. The centromere-proximal inversion O7 (yellow) is assumed to have originated along the branch of O4 . Oms became extinct as a single inversion in D. subobscura. Note that O3 is not in the path from Oa to O3+4+7 , being the inversion that generated the OST arrangement. (C) Five decades of cyclic seasonal change of O3+4+7 at Mount Pedroso, Spain. Consecutive seasonal data (dots) from the same year are connected by lines. The gray background plots the ± 2σ confidence band around the seasonal averages, and the red dot the summer-like value recorded during the spring 2011 heatwave. Included are published data from 1976 to 1981 (Fontdevila et al., 1983), 1988 to 1991 (Rodríguez-Trelles et al., 1996), 2011 to 2012 (Rodríguez-Trelles et al., 2013), and our 2015 unpublished arcsin-transformed records for late summer (0.845) and autumn (0.574). SP, spring; ES and LS, early and late summer; AU, autumn. MATERIALS AND METHODS of the flies (Zollinger, 1950), became the first established (Koske and Maynard-Smith, 1954). The chromosomal arrangements Species Karyotype and Inversion of the strain, which happened to be those most common in Central Europe, were subscripted ST (for “standard”) and from Nomenclature them new inversions were designated with numeral subindices Drosophila subobscura shows the ancestral karyotype following their order of discovery (Kunze-Mühl and Sperlich, configuration of the genus Drosophila, consisting of five large 1955). This naming system was not intended to convey polarity telocentric rods (Muller elements A-E) and one dot (Muller F) of evolutionary change. Accordingly, O3+4+7 is the arrangement (Powell, 1997). The five rods include the sex chromosome (Muller that can be interconverted with OST by the two centromere-distal A) and four autosomes of which the O chromosome (Muller E; overlapping inversions O3 and O4 (denoted by the underline homologous to chromosome arm 3R from D. melanogaster) is the joining the subscripts; Zouros et al., 1974) and the centromere- largest (∼30 Mb), comprising around 25% of the species’ nuclear proximal inversion O7 . The ancestor-descendant relationships of euchromatic genome (∼125 Mb; Karageorgiou et al., 2019). these inversions are shown in Figure 1B. An early landmark in the study of chromosomal inversion polymorphisms of D. subobscura was the development of structurally homozygous strains, as tools to identify new Drosophila Lines inversions by the location and shape of the loops formed in O chromosome conserved synteny analysis was based on data inversion heterozygotes (Zollinger, 1950; Maynard-Smith and from four whole-genome de novo assemblies, including three Maynard-Smith, 1954; Zouros et al., 1974; Loukas et al., 1979). PacBio long-read–based assemblies from D. subobscura and one The “Küsnacht” strain, named after the Swiss locality of collection Illumina short-read–based assembly from D. guanche. Of the Frontiers in Genetics | www.frontiersin.org 4 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance three D. subobscura assemblies, one was used as reference for gel electrophoresis. The genome of the Ds_7 isogenic line inversion O7 and was newly generated in this study. The other was sequenced to nominal 66-fold genome coverage using two were used as references for the standard configuration PacBio (Pacific Biosciences, Menlo Park, CA, United States) [note that the distal breakpoint of O7 maps within inversion Sequel single-molecule real-time (SMRT) technology from Oms (Karageorgiou et al., 2019), whereby is expected to exhibit a 20-kb SMRTbell template library, using Polymerase 3.0 opposite orientation in D. subobscura relative to D. guanche; chemistry and two SMRT cells. Libraries construction and Figure 1B] and were already available (Karageorgiou et al., 2019; PacBio sequencing were outsourced to Macrogen (Macrogen Inc., Bracewell et al., 2019). Also available was the assembly from Seoul, South Korea). D. guanche (Puerma et al., 2018), which was used as an outgroup. Henceforth, we will refer to these four assemblies as Ds_7, Ds_ch- Chromosome-Scale Assembly and cu, Ds_B, and Dg, respectively. Scaffolding To generate the Ds_7 assembly, we developed a line Raw PacBio reads were assembled using the Canu assembler that is isogenic for an O3+4+7 arrangement from the wild (version 1.8; Koren et al., 2017) on recommended settings for and homokaryotypic and highly homozygous for the ST read error correction, trimming and assembly, and genome size arrangements of the rest of the chromosomes (i.e., AST , JST , set at 150Mb based on previously published flow cytometry UST , EST , and O3+4+7 ). The O arrangement was first isolated data (Karageorgiou et al., 2019). These analyses were performed by crossing wild males to virgin females from the cherry-curled on a 2.80-GHz 8-CPU Intel Xeon 64-bit 32GB-RAM computer (ch-cu) recessive marker stock; they were then submitted to running Ubuntu 18.04 LTS. nine generations of backcrossing with ch-cu females and finally Chromosome-scale assembly and scaffolding followed the isogenized using the Va/Ba balancer stock (Sperlich et al., 1977). four steps outlined in Karageorgiou et al. (2019) as well as The expression of the Ba gene is highly variable. Therefore, to a fifth step, to improve genome completeness and contiguity, prevent potential errors at sorting out phenotypically O3+4+7 consisting of merging the Ds_7 assembly with a preselected set homokaryotypes, the Va/Ba stock was previously selected for of segments from the reference Ds_ch-cu assembly using one zero macrobristles on the scutum and scutellum. Crossing round of quickmerge (Chakraborty et al., 2016), as follows: first, schemes and the methods for polytene chromosome staining the CANU contigs that could be certainly anchored, ordered, and and identification are described elsewhere (Rodríguez-Trelles oriented on the nuclear chromosomes were aligned against the et al., 1996). The assayed line was stored frozen at −80◦ C Ds_ch-cu reference using NUCmer (Kurtz et al., 2004). Second, immediately upon obtention. The wild flies used to develop the the segments of Ds_ch-cu not overlapped by the CANU contigs, line were derived from our survey of the natural population each extended 10 kb outward from each of its two ends, were of Berbikiz (Spain; Lat.: 43,18949, Long.: –3,09025, Datum: extracted. Finally, separately for each chromosome, the extracted WGS84, elevation: 219 m a.s.l) conducted in July 7, 2012 Ds_ch-cu segments, together with the CANU contigs set as the (Rodríguez-Trelles et al., 2013). backbone, were fed into quickmerge. This approach was found The remaining three assemblies were derived from strains to reduce the chances of misassembly and chimerism, while homokaryotypic for all chromosomes. The Ds_ch-cu assembly making it straightforward to trace the non-backbone sequence was generated from the ch-cu strain of our laboratory (AST , in the assembly. Dot plots of the merged assembly against JST , UST , EST , and O3+4; Karageorgiou et al., 2019) and the the reference Ds_ch-cu assembly were used as a further step Ds_B assembly from an isofemale laboratory stock derived from of misassembling correction. The obtained Ds_7 assembly was a natural population from Eugene, Oregon, in 2006 (AST , JST , polished with 26 × mean coverage of 150–base-pair (bp) MP U1+2 , EST, and O3+4 ; Bracewell et al., 2019). The Dg assembly Illumina reads from the O3+4+7 isogenic line using two rounds was generated from an isofemale laboratory stock derived from a of PILON (version 1.23; Walker et al., 2014). natural population from the Canary Islands, Spain, in winter 1999 (Puerma et al., 2018); it shows the chromosome configuration Genome Annotation of the last common ancestor of the subobscura subgroup except Gene prediction and annotation of the assembled genome for chromosome E, which carries the arrangement Eg 1 (Aa , Ja , were conducted using the MAKER (version 3.01.02.-beta; Holt U1+2 , Eg 1 , and Oa ; Puerma et al., 2018; Karageorgiou et al., 2019; and Yandell, 2011; Campbell et al., 2014) annotation pipeline. Bracewell et al., 2019). Repetitive elements were identified using RepeatMasker (version 4.0.6; Smit et al., 2013/2015, at1 ) combined with three repeat High Molecular Weight Genomic DNA libraries, including (i) the Drosophila genus–specific repeat library contained in the Repbase database (release 20170127; Isolation and PacBio Whole-Genome Bao et al., 2015); (ii) a library of subobscura subgroup specific Sequencing satellites, sat290 and SGC-sat (Karageorgiou et al., 2019); and High-quality high-molecular-weight gDNA was obtained from (iii) a library of de novo identified repeats generated using 60 mg of −80◦ C frozen adult females, using a modified version RepeatModeler (version1.0.11) on the assembly masked for the of the phenol/chloroform method of Chen et al. (2010) that first two libraries. Novel long terminal repeats (LTRs), miniature yields ∼25 µg of high-quality DNA per assay, as assessed by inverted-repeat transposable elements (MITEs), tandem repeats, NanoDrop ND1000 (NanoDrop Technologies Inc., Wilmington, DE, United States) spectrophotometer and standard agarose 1 http://repeatmasker.org Frontiers in Genetics | www.frontiersin.org 5 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance and rDNA and tDNA genes were identified using LTRharvest uninverted and inverted states to isolate segments A, B, C, and (GenomeTools version1.5.10; Ellinghaus et al., 2008), MITE D, correspondingly. The remaining sequences of the uninverted Tracker (version 2.7.1; Crescente et al., 2018), Tandem Repeat state were submitted to a second round of comparative analysis Finder (TRF; version4.09; Benson, 1999), RNAmmer (version 1.2; among them, and with segments A to D to identify the Lagesen et al., 2007), and tRNAscan-SE (version 2.0; Lowe and homologies missed in the first round. As representatives of the Chan, 2016), respectively. All tools were run on default settings, uninverted state, we used Ds_ch-cu together with the previously except LTRharvest, for which we set -seed 100, -similar 90.0, and - published assemblies Ds_B and Dg, and this last one was set mintsd 5, following Hill and Betancourt (2018). The quality of the as the outgroup. annotation was controlled using the Annotation Edit Distance (AED) metric (Eilbeck et al., 2005). AED values are bounded Phylogenetic Inferences between 0 and 1. An AED value of 0 indicates perfect agreement MAFFT-based tree reconstruction of the Attacin gene family of the annotation to aligned evidence, and conversely, a value of in Drosophila was performed via maximum likelihood. Model 1 indicates no evidence support. selection and tree inference were conducted using IQ-Tree Functional annotation of MAKER-predicted proteins (Kalyaanamoorthy et al., 2017; Nguyen et al., 2015). Tree was made by BLASTP (version 2.6.0 +) searches against the searches were conducted starting from sets of 100 initial Drosophila UniProt-SwissProt manually curated datasets maximum parsimony trees using nearest neighbor interchange (Apweiler et al., 2004). Prediction of protein functional domains with default perturbation strength and a stopping rule settings. was accomplished using InterProScan (version 5.29–68.0; Jones Branch support was assessed using the ultrafast bootstrap et al., 2014) on the Pfam (Finn et al., 2016), InterPro (Finn approximation (UFboot; 1,000 replicates) (Hoang et al., 2018), et al., 2017), and Gene Ontology (Ashburner et al., 2000; The and two single-branch tests including the Shimodaira–Hasegawa- Gene Ontology Consortium, 2017) domain databases. Genome like approximate likelihood ratio test (SH-aLRT; 1,000 replicates) assembly and annotation completeness were gauged using the (Guindon et al., 2010) and the approximate Bayes parametric test Benchmarking Universal Single-Copy Orthologs (BUSCO) tool (Anisimova et al., 2011). [BUSCO, version 4 (Seppey et al., 2019)], with the latest update of the dipteran gene set (diptera_odb10), which contains 3,285 highly conserved, single-copy genes expected to be present in Non-B DNA Sequence and Transcription any dipteran genome. Factor Binding Site Scans Scans for potential non-B DNA–forming sequences considered Isolation and Characterization of the O7 the following features: inverted repeats (IRs) (capable of Breakpoints forming hairpin and/or cruciform DNA), direct/tandem repeats Suppose that +A|+B+C|+D and +A|−C−B|+D represent two (slipped/hairpin structures), mirror repeats (triplexes), alternate chromosome arrangements whose gene orders differ only by purine-pyrimidine tracts (left-handed Z-DNA), G4 motifs the orientation of the segment between A and D (with (tetraplex and G−quadruplex DNA), and A−phased repeats symbols denoting A and D, the segments upstream from (static bending). Searches were conducted online using for IRs the centromere-proximal breakpoint and downstream from Palindrome Analyzer (Brázda et al., 20163 ; accessed January 24, the centromere-distal breakpoint, respectively; vertical bars, 2020) with repeat length of 6-20 nt, spacer length ≤ 10 nt, breakpoint junctions; and plus/minus signs, orientation of the and number of mismatches ≤ 1; for tandem repeats Tandem segment relative to the uninverted sequence). We proceeded in Repeat Finder (TRF version 4.09; Benson, 19994 ; accessed Jan two steps. First, we isolated the regions containing the breakpoint 24, 2020) in basic mode; and for the remaining features nBMST junctions by chromosome conserved synteny analysis between (Cer et al., 20125 ; accessed January 24, 2020) with prefixed default the uninverted and inverted states using the Synteny Mapping settings. The propensity of IRs to adopt non-B conformation was and Analysis Program (SyMAP, version 4.2.; Soderlund et al., assessed using the difference in free energy between the DNA 2011) tool on default options, and NUCmer (see Karageorgiou sequence in the linear and cruciform structures, as implemented et al., 2019). The O7 breakpoints were identified as the loci in Palindrome Analyser (Brázda et al., 2016). of interrupted synteny whose locations and distance from each Transcription start site (TSS) prediction was conducted using other agree with the cytogenetic mapping data of the inversion the NNPP method (Reese, 20016 ). Searches for putative binding (Karageorgiou et al., 2019). Second, we localized the breakpoint sites for Relish (Rel), the heterodimer Dif/Rel, dFOXO, Dorsal junctions at base-pair resolution and performed comparative (dl), and Serpent (srp) transcription factors in the 1-kb upstream analyses of their flanking sequences using the progressive guide region of the Attacin predicted TSSs were performed using tree-based MAFFT algorithm (version 72 ) with the accuracy- the FIMO tool (Grant et al., 2011) from the MEME suite oriented method “L-INS-i” (Katoh et al., 2019). Each of the (Bailey et al., 2015). For Rel and Dif/Rel, and for dFOXO, regions +A|+B and +C|+D from the uninverted state was aligned separately, first with +A|−C and then with −B|+D from the 3 inverted state. From each of the four resulting alignments, we http://bioinformatics.ibp.cz:9999/#/en/palindrome 4 used the regions showing positional homology between the https://tandem.bu.edu/trf/trf.html 5 https://nonb-abcc.ncifcrf.gov/apps/nBMST/default 2 6 http://mafft.cbrc.jp/alignment/software/ https://www.fruitfly.org/seq_tools/promoter.html Frontiers in Genetics | www.frontiersin.org 6 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance we used the FootprintDB database (Sebastián and Contreras- Moreira, 20147 ) Drosophila melanogaster Major Position Matrix Motifs (DMMPMM) identified, respectively, by Senger et al. (2004) and Weirauch et al. (2014). For dl and srp, we used the REDfly database (version 5.5.3; Rivera et al., 20198 ) improved iDMMPMM motifs developed by Kulakovskiy and Makeev (2009). Searches were performed using a p value cutoff of 10−3 . RESULTS Chromosome-Scale Assembly and Annotation of Chromosome Arrangement O3+4+7 The PacBio Sequel sequencing of the O3+4+7 isogenic line genome generated 2,457,493 reads, with mean and longest lengths of 11,257 bp and 117,750 bp, respectively. Canu correction and trimming retained a 42-fold genome coverage for the assembly. Of the 385 Canu-generated contigs, the 14 that could be confidently anchored, ordered, and oriented covered the complete reference genome, with an added length of 126.770 Mb and N50 of 10.587 Mb. Quickmerge of those 14 CANU contigs resulted in six chromosome-scale scaffolds, one per each of the major D. subobscura chromosomes (Table 1). Of note, chromosome O was built from two contigs only, with the centromere-proximal contig (tig00026085; 29.679 Mb) spanning almost all the chromosome length (96.9%) (Figure 2A). The Ds_7 assembly contained 13,459 MAKER-annotated genes, nearly all with well-supported predictions (AED50 = 99.3%). Only 2.6% (87) of the BUSCO genes were missing, indicating that the assembly is almost complete. The O chromosome contained FIGURE 2 | Chromosome conserved synteny analysis of the O7 breakpoints. 3,220 (23.9%) of the annotations of the assembly. (A) The long read-based O3+4+7 chromosome-scale scaffold (s001). The vertical dotted line near the telomere indicates the location of the stitch Identification of Inversion O7 Using between the two Canu tigs. (B) SyMAP comparative synteny analysis showing the inversions found along the path from the ancestral arrangement of the Chromosome Conserved Synteny subobscura subgroup (Oa ) to O3+4+7 (see Figure 1B). In addition to O7 , the Analysis two chromosomes also differ by inversions O4 and Oms . (C) SyMAP direct comparison of the Oa and O3+4+7 chromosome arrangements. Bands The structural transition between the O chromosomes of the connecting the chromosomes denote uninverted (pink) and inverted (green) Ds_7 and Ds_ch-cu assemblies called for one large megabase- synteny blocks. Labeled ticks on chromosomes indicate proximal (p) and sized inversion (Figures 2B,C), whose breakpoints located distal (d) inversion breakpoints. The remaining symbols are as in Figure 1. 7 http://floresta.eead.csic.es/footprintdb 8 http://redfly.ccr.buffalo.edu/ cytologically precisely as it would be expected if they were from O7 . Relative to the nearest of the available 140 cytologically mapped markers of the O chromosome (see Karageorgiou TABLE 1 | Ds_7 assembly summary statistics (Muller elements are given in parenthesis, and lengths are given in megabases of sequence). et al., 2019), the proximal breakpoint was located 44.5 kb downstream from Sb (Dmel\CG4316) and 117.4-kb upstream Component Length Scaffolds Canu Largest Canu Gene from microsatellite dsub02, and the distal breakpoint 111.8 kb contigs contig models downstream from rdx (Dmel\CG12537) and 29.3kb upstream Nuclear 126.770 6 14 29.679 13,459 from Abi (Dmel\CG9749). Sb and dsub02 have been respectively genome mapped to subsections 77B (Dolgova, 2013) and 77C (Santos Dot (F) 1.412 1 1 1.412 96 et al., 2010), and rdx and Abi to subsection 85E (Dolgova, 2013; A (A) 22.941 1 2 17.229 2,323 Pegueroles et al., 2013) of the Kunze-Mühl and Müller (1958) J (D) 25.018 1 3 10.587 2,652 standard cytological map. Other than O7 , no D. subobscura U (B) 26.010 1 3 13.133 2,561 inversion maps to those positions. E (C) 20.783 1 3 9.524 2,607 Comparative analysis of the genes annotated in the regions O (E) 30.629 1 2 29.679 3,220 immediately flanking the breakpoints in Ds_7, Ds_ch-cu and Frontiers in Genetics | www.frontiersin.org 7 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance FIGURE 3 | Sequence-annotated breakpoint regions in O7 and the standard uninverted state. (A) Gene scale (10-kb scale bar) depiction of the proximal and distal breakpoints in the uninverted state (+A|+B and +C|+D, respectively) versus O7 (+A|−C and −B|+D). In green are segments A and B, and in sepia C and D. Vertical broken lines indicate break junctions, and arrow boxes the size and direction of the genes labeled vertically using the names of the corresponding D. melanogaster orthologs. The coordinates of O7 in the assembled chromosome are given in parentheses. (B) Zoom-in (1-kb scale bar) on the regions immediately flanking the break junctions, with O7 oriented backward (i.e., +B+C, instead of −C−B) to better track the differences with the uninverted state. Arrow boxes indicate the size and direction of the sequence elements discussed in the text. Gray boxes (exons) linked by polygonal lines (introns) represent the two AttA2 paralogs oriented in the direction of transcription. In the distal breakpoint, the two alternative haplotypes of the uninverted states, namely, that from ch-cu and that from Dg and B, are shown. Note the reversal of the spacer in ch-cu versus Dg and B, and the mirror halves flanking dd7 in O7 . (C) O7 represented as in (B), but in its actual orientation (i.e., -C-B). Ds_B with those in the outgroup Dg (Figure 3A) corroborated in Supplementary Figures 1, 2. Figure 3B provides a schematic that Ds_ch-cu and Ds_B carried the uninverted state, whereas representation of the +A|+B region based on the alignment Ds_7 carried the inverted state. The assembled O7 has a size of of Supplementary Figure 3. In the case of O7 , the region 9,936,431 bp, totaling 32.4% of the chromosome (30,629,152 bp). was reconstructed using the reverse complement of segment - It has a GC content (43.8%) below that of the O chromosome B. The breakpoint junction is located within a 2,445-bp-long (44.9%) since it is located in the chromosome centromere- sequence stretch present only in the inverted state. The site of proximal half, which is relatively AT-rich (Karageorgiou et al., the insertion is flanked by multiple indels, which suggests that 2019). O7 was predicted to have 1,028 protein-coding genes, the insertion occurred in a region of prior sequence instability. or 31.9% of the gene models of the O chromosome, in close Of the insertion length, 2,317 bp are on the + A segment and agreement with its percent of chromosome length. 128 bp on the + B segment. The insertion begins with a 153-bp- long direct repetition (R1-2) of the upstream flank. Proceeding downstream from this repeat, there are two inverted duplications Nature and Properties of the DNA named d1 and d2, each with copies A and B, with d1 shorter Sequences Surrounding O7 Breakpoint (59 bp long each of d1A and d1B) than d2 (534 and 540 bp for Junctions copies d2A and d2B, respectively). The two A copies (i.e., d1A and Proximal Breakpoint of O7 d2A) are separated from the two B copies (i.e., d1B and d2B) by The alignments used for isolation of the breakpoint junctions and an intervening sequence of 1,374 bp. The junction between + A their corresponding flanking regions A, B, C, and D are shown and + B is precisely located between d1B and d2B. d2B extends Frontiers in Genetics | www.frontiersin.org 8 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance FIGURE 4 | Isochromatid-chromosome mixed staggered model of O7 formation. (A) Start from an individual heterozygous (or homozygous) for the insertion at the proximal breakpoint, and heterozygous ch-cu–type (above)/B-type (below) for the orientation of the spacer at the distal breakpoint. (B) Occurrence of two pairs of isochromatidal staggered SSBs at the proximal breakpoint (open squares), and two staggered chromosomal blunt-ended DSBs at the distal breakpoint (open rectangles), with demarcation of breakpoints flanking segments +A|+B and +C|+D, respectively. (C) Emergence of broken ends with single stranded overhangs (+ d1A, + d1B and + d2A, + d2B) at the proximal breakpoint. (D) Microinversion formation with fill-in of the overhangs, resulting in the terminal inverted duplications + d1A and -d1B, and the reversed orientation of segment -d2A. And O7 formation by reversal of the +B+C segment with fill-in of the -d2B overhang, and distal breakpoint repair via rejoining in trans with the homologous chromosome. 409 bp downward from the downstream end of the insertion into Relative to the predicted nearest gene TSSs, the events the region of resumed homology between O7 and the uninverted took place in an intergenic region. Specifically, the upstream- state, indicating the orientation of the parent copy. most SSB occurred 1,364 bp downstream from Akt1 (CG4006; The above pattern of sequence copy number, order, and serine/threonine–protein kinase B) gene, and the downstream- orientation most parsimoniously indicates that the proximal most one 2,047 nt upstream from Mhcl (CG31045; myosin breakpoint of O7 was formed on an insertion region that heavy chain-like) gene (Figure 3A). From our repeat annotation experienced two pairs of staggered SSBs, which resulted in pipeline, the region around the breakages is a composite of two DSBs (Figure 4; but see section “DISCUSSION” for repetitive sequences [16 in total, ranging in length from 21 bp alternative models). The upstream-most DSB generated the of a (TTG)n simple-repeat to 532 bp of satellite rnd-4_family- proximal breakpoint of a 2,026-bp-long microinversion and the 179], interspersed with traces of transposable elements [84 bp downstream-most DSB generated a junction flanked upstream from an LTR and 72 bp from a long interspersed nuclear element by the distal end of the microinversion and downstream by (LINE)]. Overall, no evidence of open reading frames and/or the proximal end of O7 . Accordingly, duplications d1 and d2 specific motifs could be found pointing to the observed breakages would, respectively, represent the filled-in staggered SSB-induced as directly caused by the insertion/excision of other sequences. terminal single-stranded overhangs of the microinversion and The role of non-B DNA as source of DSBs is well-established. inversion O7 . Figure 3C shows O7 ’s segments A and B such as Generally, DSBs are expected to colocalize with their causal they are found in the inversion. That d2A and d2B show direct non-B DNA motifs (e.g., Kolb et al., 2009; Lu et al., 2015). instead of reverse relative orientation as it would be expected if We used this prediction to investigate whether the local DNA paired–staggered SSBs generate inversions with inverted repeated conformational environment of the ancestral sequence could ends (Ranz et al., 2007) would be explained by the reversal in the have acted as trigger or mediator of the complex rearrangement orientation of d2A as a result of the microinversion. of the proximal breakpoint region. We proceeded in two steps: Frontiers in Genetics | www.frontiersin.org 9 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance of 8.7 copies of the consensus heptanucleotide AATAAAT, and one mirror repeat of two 11 nt-long repeats separated by a 30-nt spacer, indicating that the proximal breakpoint of O7 occurred on an unstable sequence with potential for adopting multiple alternative non-B DNA conformations. Distal Breakpoint of Inversion O7 Figure 3B provides a schematic representation of the +C|+D region based on the alignment of Supplementary Figure 6. In the case of O7 , the region was reconstructed using the reverse complement of segment -C. From up to downstream, the breakpoint junction is located within a 450-aligned-sites-long gap-rich spacer region, spanning between two highly identical long IRs, IR1 and IR2, of 1,117 and 1,112 sites of alignment length, respectively. There is no evidence of duplicated sequence in Ds_7 relative to the other assemblies, indicating that the DSB either was a clean cut or did not involve significantly staggered FIGURE 5 | Distribution of inverted repeats (IRs) with potential to adopt non-B SSBs. On the other hand, the spacer of Ds_7 was the shortest DNA hairpin/cruciform structures across a 10-kb window centered at the region of the O7 proximal breakpoint, as obtained using Palindrome Analyzer (250 nt) of all four lines (407, 317, and 343 nt for Ds_ch-cu, (Brázda et al., 2016). The region of interest is shown at scale above the plot. Ds_B, and Dg, respectively) because of a single deletion located The purple and red boxes represent respectively the parental copies of d1 and precisely at the center of the region (hereon called dd7, for d2 prior to their duplication as a result of the two pairs of staggered SSBs distal deletion of O7 ). A closer look at the pattern of pairwise (black arrowheads). The highest concentration of IRs occurs around the junction of the O7 breakpoint (+A|+B). sequence similarities along the spacer revealed two findings: (i) dd7 split the Ds_7 spacer in two mirror halves. For the upstream half, Ds_7 is almost identical (96.8%) to Ds_ch-cu while bearing first, we reconstructed the region of the rearrangement before no detectable homology to Ds_B, whereas for the downstream the breakages. It should be recalled that most of the rearranged half, Ds_7 is almost identical (97,6%) to Ds_B while bearing no sequence is embedded in an insertion that is absent in the detectable homology to Ds_ch-cu; and (ii) the spacer of Ds_ch- ancestral non-rearranged state. Therefore, we reconstructed the cu is almost identical (95.4%; excluding indels) to that of Ds_B prebreakages state by undoing the hypothetical rearrangement but in reversed orientation. The reversal occurred in Ds_ch-cu, steps that generated the present sequence state. Specifically, we because in Ds_B the spacer is oriented as in the outgroup Dg. reversed the orientation of the microinversion (Supplementary Altogether, the above observations can be understood as Figure 4) and deleted one copy of each DSB-induced duplication follows (Figure 4). Prior to the origination of the distal (Supplementary Figure 5). The resulting sequence had the breakpoint of O7 , a carrier of an uninverted chromosome of form: + d1, (+ 1,374 bp), |, + d2 (Figure 4B). Which copy of B-type experienced a reversal of the spacer region between the each of the two duplications to eliminate was inconsequential, IRs, giving rise to the uninverted chromosome of ch-cu–type. because they are nearly identical to each other in the two cases Later on, a homokaryotype for the uninverted chromosome (98.3% and 95.6%, for the identities between copies A and B that was heterozygous for the microinversion of the spacer of dup1 and dup2, respectively). Furthermore, the observed underwent at least two DSBs, one in each of two homologous high level of identity (97.3%) between d2 and its homologous non-sister chromatids, such that the DSB in the ch-cu–type region in Ds_ch-cu and Ds_B suggested that the rearrangement is chromatid occurred immediately before the first site of the recent enough to allow assuming that the original conformational dd7 and that in the B-type chromatid immediately after the sequence features that could have mediated it are still observable. last site of the dd7. Finally, the reversed + B end generated After establishing the prebreakage sequence, we next looked for by the proximal staggered DSB in the ch-cu–type chromatid sequences with the potential to form non-B DNA structures along illegitimately joined with the + D end generated by the distal DSB a 10-kb window centered on it. in its homologous non-sister B-type chromatid, which resulted in Figure 5 shows the distribution of the number of IRs capable a recombinant chromosome carrying the inversion O7 with the of forming hairpin and cruciform structures along the target exact observed dd7 deletion. sequence. The highest density occurs immediately around the Like the proximal breakpoint, the distal breakpoint occurred junction between the microinversion and inversion O7 . In in an intergenic region yet at comparatively much shorter particular, the breakpoint is located within a ∼150-nt-long distance (∼390 bp) to the nearest genes. Specifically, the breakage stretch of AT-rich sequence [simple repeat (ATTT)n, from our separated two copies of an Attacin gene (CG10146; AttA) located genome annotation pipeline] containing 15 IRs, of which one opposite to each other on each of the two arms of the long IR. Our located 68 nt downstream the breakpoint junction ranked in repeat annotation pipeline did not identify repetitive sequences in the top 5% with highest likelihood of intrastrand annealing to the vicinity of the distal breakpoint in Ds_ch-cu or Ds_B. form a hairpin (AATTTT AAAATT; 1GS – 1GL = 2.64). In We searched the region of the spacer for potential non- addition, embedded in the IR cluster, there is one tandem repeat B DNA–forming sequences in the vicinity of the breakpoint Frontiers in Genetics | www.frontiersin.org 10 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance FIGURE 6 | Pre-O7 history of instability of the distal breakpoint. (A) In the most recent common ancestor of the Drosophila genus, AttD was the only Attacin gene present in Muller element E. (B) Later, the ancestor of the obscura group lost AttC, and underwent DNA-based interchromosomal transposition of AttA (or a close paralog; see Supplementary Figure S7) from Muller C to Muller E, followed by DNA-based intrachromosomal transposition within Muller E, giving rise to AttA2 and AttA3 (whether simultaneously or sequentially and, if the latter, which was first is unknown). (C) Before the split of the subobscura subgroup, AttA2 was duplicated, giving rise to the inverted duplicates AttA2a (parent copy) and AttA2b (daughter) separated by a short central spacer. (D) In D. subobscura, the central spacer underwent a reversal, generating a microinversion polymorphism with segregating states B-type (ancestral) and ch-cu–type (derived). Genes are represented as solid black boxes (exons) linked by polygonal lines (introns), and oriented in the direction of transcription. The central spacer is represented as a box colored in three shades of gray pointing in the direction of its orientation. junctions in Ds_ch-cu and Ds_B. In both cases, we found that run closely associated with a highly dynamic evolution of the the IR with the highest propensity to form a hairpin was a perfect Attacin immunity gene family in the obscura species group. 14-bp-long palindromic sequence located next to the breakpoint This conclusion is based on phylogenetic analysis of the Attacin junctions (ATGAACT AGTTCAT; 1GS – 1GL = = 2.05; family in Drosophila (Supplementary Figure 7) using synteny to located 13 and 2 bp upstream and downstream the junction in distinguish orthologous from paralogous copies (Supplementary Ds_ch-cu and Ds_B, respectively). Apart from IRs, we did not Table 1). The results are summarized in Figures 6A–D. The most detect additional potential non-B DNA sequences around the recent common ancestor of the Drosophila genus (Figure 6A) distal breakpoint. carried three copies of the gene with relationships [(A,C),D], All nucleotides in the +A|−C region of Ds_7 could be of which the more distant D was located in Muller element unambiguously ascribed to segment A or C. However, in E, and the closer to each other A and C in Muller element the −B|+D region −B and +D are separated by 21 extra C. After it split from the melanogaster group (Figure 6B), inserted nucleotides (i.e., GAGCACTCTCCACAGCAAAGT). the branch leading to the obscura group lost copy C and We decided to ascribe this sequence to the distal breakpoint underwent an interchromosomal transposition of copy A from junction, because it contains an 8-bp substring (underlined) that Muller element C to E. The daughter copy then underwent resembles the beginning of the +D end (CATCAAAG), and hence another, in this case intrachromosomal, transposition, which it likely represents filler DNA generated by a microhomology- originated two new Attacin copies that we called AttA2 and templated repair mechanism. AttA3, with AttA2 located between foxo and Npc2b, and AttA3 located ∼300 kb downstream from AttA2, between Cul5 and Pre-inversion Record of Rearrangement Sirt7. The two transpositions were genome-based duplications of O7 Breakpoints rather than retroposition events, because the new copies Previously, it was shown that the proximal breakage of O7 was conserved the intron position of their parental gene. Before preceded by an insertion. Likewise, the region of the distal the split of the subobscura subgroup (Figure 6C), copy AttA2 breakage had a pre-inversion history of rearrangement, which underwent an inverted duplication that generated the two closely Frontiers in Genetics | www.frontiersin.org 11 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance spaced copies AttA2b and AttA2a in head-to-head orientation, and transcribed in opposite directions. In D. subobscura (Figure 6D), the spacer between the IRs experienced a reversal of orientation generating the microinversion polymorphism of the distal breakpoint. Subsequently, a heterozygote for the microinversion underwent distal DSBs that allowed the formation of the recombinant O7 inversion via ectopic repair of non-sister chromatids. Potentially Functional Effects of the O7 Mutation The distal break of O7 disrupted concerted evolution between two subobscura subgroup-specific AttA2 duplicates. This conclusion is based on the previous section’s results, together with the phylonetwork of coding sequences shown in Figure 7. FIGURE 7 | Phylogenetic network of AttA2 orthologous and paralogous Accordingly, right after the duplication of AttA2, the two sequences from the inverted (O7 ) and uninverted (B and ch-cu) states in paralogs began to evolve in concert, converting each other D. guanche and D. subobscura. The duplication is older than the species, but to generate their present characteristic phylogenetic pattern of the paralogs cluster within species owing to concerted evolution. Depicted greater resemblance between paralogs from the same species (i.e., thicker is the branch leading to the copy relocated by O7 (i.e., AttA2b), which is longer than that leading to the copy that remained in place (AttA2a) due to D. guanche and D. subobscura) than between orthologs from an acceleration of the synonymous substitution rate, likely as a result of having different species (e.g., Puig-Giribets et al., 2019). At one end of escaped concerted evolution. The split network was constructed using the the resemblance, it is the ch-cu strain, whose two AttA2 copies are NeighborNet method as implemented in SPLITSTREE version 4.14.5 (Huson identical to each other, and at the other end O7 , where the copy and Bryant, 2006), on the JC69 + I (% of invariable sites 81.6) best-fit model relocated by the inversion evolved significantly faster than the one distances obtained using the DIVEIN web server (https://indra.mullins.microbiol.washington.edu/DIVEIN/) (Deng et al., 2010). that remained in place, owing exclusively to an acceleration of Sets of parallel edges represent conflicting topological signals. the synonymous substitution rate [P < 0.05; Tajima’s relative rate test (Tajima, 1993) using either of the remaining six sequences as outgroup], as the two copies are identical at the amino acid level. the AttA2 genes of the inverted chromosome had TFBSs for The acceleration took place in the direction of a slight decrease the metabolic factor dFOXO. Furthermore, the dFOXO TFBSs in codon bias in the relocated copy (Nc = 51.2 vs. 50.7, for were all contributed by the newly attached sequence. The fact the comparison AttA2b vs. AttA2a, respectively; where Nc is the that the AttA2 genes were conserved at the amino acid level in improved effective number of codons index; Sun et al., 2013). The D. subobscura, together with the observed qualitative difference increased synonymous rate can be understood, in part because in predicted cis-acting sequence between uninverted and inverted the inversion released the two Attacin copies from evolving chromosomes, suggests that the inversion O7 brought the AttA2 in concert; and in part assuming that the expression of the genes under the influence of the IIS metabolic pathway. paralogs shifted as a result of changes in regulatory environment In addition to the Attacin immunity genes, the breakpoint associated with their relocation. regions include Akt1 and foxo, two interacting core components Considering the short spacing between the two AttA2 paralogs of the IIS metabolic pathway identified by other studies as in the uninverted chromosome (∼390 bp), it appeared likely candidate for climate adaptation (Fabian et al., 2012; Paaby et al., that the inversion would have detached them from part of 2014; Kapun et al., 2016; Durmaz et al., 2019). The roles of these their promoters, binding them to new potentially cis-acting genes and the potential impact of O7 on them are dealt with in elements. To assess this possibility, we searched 1 kb upstream the Discussion. of the predicted TSS of each gene for putative transcription factor binding sites (TFBSs) for five transcription factors (TFs), including the nuclear factor κB factors dorsal (dl), dorsal- immunity related factor Dif and Relish (Rel), the GATA factor DISCUSSION Serpent (srp), and the forkhead factor dFOXO. The first four TFs are under control of the Toll and immune deficiency (IMD) Molecular Mechanism of O7 Formation immunity pathways and regulate Attacin inducible expression in O7 Is a Complex Multibreak Inversion Formed via response to bacterial infection (Senger et al., 2004). dFOXO TF is Rejoining in trans With the Two Homologous controlled by the insulin/insulin-like growth factor signaling (IIS) Chromosomes metabolic pathway and regulates constitutive Attacin expression Sequence data on inversion formation in Drosophila have been in non-infected flies suffering from energy shortage or stress interpreted in terms of two major mechanisms with associated (Becker et al., 2010). The results are shown in Figure 8. The distinctive footprints. The first mechanism is intrachromatidal AttA2 genes had predicted TFBSs for the immunity related factors NAHR between inversely oriented repeats. This mechanism in both uninverted and inverted chromosome states, but only generates inversions with duplications at their ends in both the Frontiers in Genetics | www.frontiersin.org 12 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance FIGURE 8 | New dFOXO binding sites. Schematic representation of up to 1-kb sequence up and downstream predicted TSSs of AttA and the inverted duplicates of AttA2 (represented as in Figure 6) in O7 and the uninverted (B and ch-cu) states, including also the nearest flanking genes. Colored lines connecting genes designate the following: orange, region of the inverted repeats; dark and light blue, first and second halves of the spacer of the inverted repeats, respectively, oriented as the arrowheads; green and brown, the novel sequences to which the AttA2 copies became reattached by O7 , with corresponding breakpoints (+A|−C and −B|+D) indicated. The inverted repeats of B and ch-cu are folded over each other. Putative TFBs are symbolized: gray arrowheads and circles (palindromic sites), for Dorsal, Relish, and Diff/Relish; blue arrowheads for Serpent, and red boxes with an asterisk for dFOXO, respectively. Only AttA and the two Att2A copies of O7 have TFBs for dFOXO. inverted and uninverted states (Cáceres et al., 1999), which is not addition to the previously proposed intrasister and intersister the case of O7 . chromatidal exchanges. The second mechanism is chromosomal breakage and ectopic repair via NHEJ. This mechanism either does not The Breaks of the O7 Inversion Were Likely Induced generate duplications or generates them but at the ends of by Non-B DNA Secondary Structures the inverted state only. These two types of NHEJ footprints Inversion O7 provides, to our knowledge, the first compelling have been explained in terms of two alternative modes evidence for a role of non-B DNA in inversion formation in of breakage: cut-and-paste via clean DSBs that generate Drosophila. Previous studies had reported the presence of AT- blunt ends and staggered on the same (isochromatidal) or rich sequences around the breakpoints of some fixed (Cirera different (chromatidal) sister chromatids (see Introduction). et al., 1995; Richards et al., 2005) and polymorphic (Prazeres In the case of O7 , it is not a cut-and-paste inversion, but da Costa et al., 2009) inversions. In no instance, however, were neither is it a typical staggered breaks inversion. Thus, particular sequences susceptible to adopt secondary structures while the inversion proximal breakpoint could be either identified. In the case of O7 , the proximal break junction occurred isochromatidal (Figure 4) or chromatidal (Figure 9), the just within a palindromic AT-rich repeat capable of adopting distal breakpoint has to involve the two homologous hairpin/cruciform, slipped and triplex DNA conformations. chromosomes (Figures 4, 9). This latter pattern could be Likewise, the distal junctions are located next to perfect 14-bp- deduced because of the chanceful circumstance that our long hairpin/cruciform-forming palindromes. two representatives of the uninverted state (i.e., Ds_ch- The role of non-B DNA-forming sequences in causing cu and Ds_B) segregated for the microinversion of the genome instability is well-established (Wang and Vasquez, 2006; spacer between the IRs flanking the distal breakpoint. Lobachev et al., 2007; Aguilera and Gómez-González, 2008; Zhao Alternatively, the distal breakage could have occurred in et al., 2010). The shift from B to non-B DNA conformation occurs a recombinant between chromosome types ch-cu and B. while DNA is in single-stranded form, e.g., behind replication This, however, appears unlikely because crossover within forks, between Okazaki fragments, or in actively transcribed microinversions should be extremely rare (Greig, 2007). genes (Voineagu et al., 2008). Non B-DNA structures induce Our conclusion agrees with a study of the genealogical DSBs through, e.g., stalling replication and transcription (Mani relationships between inversions of the E chromosome and Chinnaiyan, 2010; Kaushal and Freudenreich, 2019). There in D. subobscura, which proposed that E9 arose in a are no specific predictions as to the type, number, and location heterokaryotype EST /E1+2 to accommodate a conflict between of the DSBs generated by any given structure in any particular molecular and cytological data (Orengo et al., 2019). This situation. Still, a single structure can induce multiple DSBs across and our results indicate that NHEJ inversions form through hundreds of base pairs around it (Wang et al., 2006; McKinney mechanisms that can incorporate information from the et al., 2020), and stalled replication forks can accumulate up to two homologous chromosomes (chromosome model), in 3 kb of single-stranded DNA (Sogo et al., 2002; Lopes et al., Frontiers in Genetics | www.frontiersin.org 13 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance FIGURE 9 | Chromosome model of O7 formation. (A) Start from an individual homozygous for the insertion at the proximal breakpoint, and as in Figure 4 at the distal breakpoint. (B) Occurrence of two pairs of staggered chromosomal blunt-ended DSBs (open rectangles) at the proximal breakpoint, and as in Figure 3 at the distal breakpoint. (C,D) Microinversion formation, and formation of O7 via rejoining in trans with the homologous chromosome as indicated. The model results in an order of the duplications (i.e., + d1A, -d2A, -d1B, -d2B) identical to that resulting from the model of Figure 4. 2006). In the case of O7 , this length is well over the size of the staggering models could in fact be footprints of a replication- overhangs that would be generated by an isochromatid model based mode of repair (reviewed in Kramara et al., 2018; Scully of the proximal breakpoint (58 and 534 nt; see Figure 3 and et al., 2019), including (i) presence of non-B DNA-forming Supplementary Figure 3). sequences just in, or adjacent to breakpoint junctions (see below); (ii) spatial proximity of the breakpoint regions in the nucleus, The Inverted Duplications at the O7 Breakpoints as evinced by the fact that the genes flanking the junctions are Could Be Footprints of Repair Instead of Staggered closely related functionally (Farré et al., 2015; but see Sunder Breakage and Wilson, 2019); and (iii) multiple breaks concentrated in All the aforementioned inverted duplication-generating NHEJ a short sequence segment. A fourth feature, namely, presence models are predicated upon the role of DNA breakage (Ranz of microhomology at the distal breakpoint junction, would be et al., 2007). However, the inverted duplications at the ends also consistent with a recombination-based mechanism such as of O7 could also be explained as a result exclusively of repair, MMEJ. Overall, these features suggest that O7 arose as result with no need for invoking staggering of the breaks. DNA repair of a non-B DNA-induced replication impairment, affecting at has emerged as a key factor capable of generating extremely least its proximal breakpoint. It is known that this type of complex breakpoint sequence rearrangements (reviewed in events can trigger BIR and MMBIR repair (Sakofsky et al., Scully et al., 2019). The spectrum of known error-prone repair 2015). Of the two pathways, the second pathway has yet to be mechanisms can be grossly classified as recombination-based, identified in Drosophila (Alexander et al., 2016; Bhandari et al., such as microhomology-mediated end-joining (MMEJ), and 2019). A possible scenario is detailed in Figure 10: first, non- replication-based, such as break-induced replication (BIR) and B DNA-induced stalling of a replication fork at the proximal microhomology-mediated BIR (MMBIR) (Lee et al., 2007; Zhang breakpoint of a ch-cu–type chromosome led to two DSBs et al., 2009; Hastings et al., 2009). Here, the term microhomology generating three fragments. Second, the centromere-proximal is used to mean a short tract (∼1 – 25 bp) of chance similarity, fragment engaged in a BIR event using the homologous region rather than common descent. In the case of O7 , three features of a B-type chromosome. Third, a second fork stalling triggered suggest that what appear to be footprints of breakage by the a switch from BIR to MMBIR with template switching to Frontiers in Genetics | www.frontiersin.org 14 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance FIGURE 10 | Chromosome and BIR/MMBIR repair model of O7 formation. (A) Start as in Figure 9. (B) Occurrence of one pair of blunt-ended DSBs (open rectangles) at the proximal breakpoint of the ch-cu type chromosome, and as in Figure 4 at the distal breakpoint. (C,D) Step 1: 50 to 30 resection generating a 30 single stranded + d1A end. Step 2: beginning of a BIR event via strand invasion into the homologous region of the B-type chromosome. Step 3: switch from BIR to MMBIR, with forward template switching to the distal end of + d2A and backward copying. Step 4: MMEJ to the distal break-end of O7 from the original ch-cu–type chromosome. The distal breakpoint repaired as in Figure 9. The model results in an order of the duplications (i.e., + d1A, -d2A, -d1B, -d2B) identical to that resulting from the model of Figures 4, 9. a downstream microhomology. Copying backward from the A gene that are highly conserved. It includes at least three new template resulted in the rearrangement of the proximal rearrangements that occurred in the lineage of D. subobscura breakpoint, including the inverted duplication of the O7 end after its separation from that of the melanogaster group (see (e.g., Lee et al., 2007; Smith et al., 2007; Carvalho et al., 2015; section “RESULTS”; Figure 6), namely, (i) insertion of AttA2 Tremblay-Belzile et al., 2015). Finally, the event was terminated between the foxo and Npc2b genes; (ii) emergence of the IR by by an MMEJ to the distal break-end of O7 from the original inverted duplication of the parental AttA2 (Figure 6B), which ch-cu-type chromosome (e.g., Scully et al., 2019). could have occurred through an event of forward template switching and backward copying by the DNA polymerase (Smith The O7 Breakpoints Carry a Pre-inversion Record of et al., 2007; Lee et al., 2007), as discussed above; and (iii) Fragility emergence of the ch-cu–type chromosome via inversion of the The breakpoint sequences of O7 had a record of instability spacer between the IRs in a B-type chromosome (Figure 6D), prior to the origin of the inversion, as evinced by the fact that which could be explained as an outcome of a stem-loop they are located within sequences inserted from elsewhere in formation by the IR, followed by resolution of the strand- the genome. This suggests that the regions that gained those exchange junctions between the IR arms (see Figure 4 in Leigh insertions were relatively exposed in the nucleus (reviewed in Brown and Ish-Horowicz, 1981; Figure 3 in Kolb et al., 2009 and Farré et al., 2015). In the case of the proximal breakpoint, that Zhao et al., 2010). could be associated with high levels of transcriptional activity at The pre-O7 insertion in the proximal breakpoint is specific the broadly expressed Akt1 gene (Andjelković et al., 1995; Slade to D. subobscura and is therefore much more recent than that and Staveley, 2016). of AttA2 in the distal breakpoint. Preliminary analyses indicate That the O7 junctions arose in fragile regions, beyond that it is internally rearranged relative to other paralogous copies, the proximate effects of their associated non-B DNA (see supporting that it carries recombinogenic potential. The origin above), may be most apparent from the pre-inversion record and evolution of this inserted sequence, as well as its possible of recurrent rearrangement of the IR at the distal breakpoint implication in the formation of other D. subobscura inversions, (Figure 6). This record is particularly amenable to reconstruction warrant further investigation (CK, RT, and FR-T; manuscript because the IR largely consists of two copies of the Attacin in preparation). Frontiers in Genetics | www.frontiersin.org 15 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance O7 Breakpoints Potentially Functional midsummer, coinciding with the growth season, it seems more Effects likely that whatever the effect of the inversion mutation on Akt1 and/or foxo, it occurred in the direction of an enhanced O7 Relocated foxo in Tight Linkage Association With basal IIS versus dFOXO activity relative to the OST ancestral Its Antagonistic Regulatory Partner of the IIS state. This would raise the question of why the O7 frequencies Metabolic Pathway Akt1 decrease (and those of OST increase) every year from late O7 changed foxo from being megabases (∼10 Mb) away summer to winter. from Akt1 to being tightly linked to it, with only the short AttA2b gene sandwiched between them. Akt1 and foxo are O7 Disrupted the Concerted Evolution of Two AttA2 functionally conserved genes, which, in Drosophila, encode Immunity Genes and Reattached Them to Putative the serine/threonine–protein kinase B AKT/PKB, and the dFOXO Metabolic Enhancers forkhead-box DNA-binding domain-containing TF dFOXO, The immune function is highly energy demanding in terms respectively. The two genes are key antagonistic regulators of both maintenance and, especially, rapid deployment upon of the IIS pathway (Teleman, 2010; Slade and Staveley, infection (reviewed in Dolezal et al., 2019). Therefore, within 2016), a major trigger of shifts in anabolic versus catabolic a limited energy budget, a trade-off is expected between cellular activity in response to nutritional status (de Jong reproduction and immunity (Schwenke et al., 2016). The and Bochdanovits, 2003) and multiple other cues (Regan Drosophila innate immune response consists of a cellular and et al., 2020). In abundant nutrient conditions, AKT/PKB a humoral component. The humoral component involves the inactivates dFOXO, thus shifting food energy allocation toward production of antimicrobial peptides, among which Attacins are reproduction and growth (the IIS pathway). Conversely, active against gram-negative bacteria (Hanson and Lemaitre, scarce nutrient conditions prevent AKT/PKB from inactivating 2020). The two main modes of Attacin production, including dFOXO, which redirects metabolism toward mobilization the induced (by a factor of even > 100) upon infection mode, of energy stores for somatic maintenance (FOXO pathway). and the basal in absence-of-infection mode link immunity Laboratory research using large effect mutants has shown with the Akt1/foxo IIS metabolic signaling pathway (Becker that the IIS/FOXO pathway is extensively pleiotropic, with et al., 2010; Dolezal et al., 2019). The inducible mode is major evolutionary conserved effects on fitness-related life- regulated primarily by the immunodeficiency Imd signaling history traits, including growth, size, reproduction, lifespan, pathway and to a lesser extent by the Toll signaling pathway. and stress resistance (reviewed in Flatt and Partridge, 2018). The two signaling pathways have the same effect of activating Research from the field found IIS loci to harbor substantial dFOXO, thus mobilizing resources toward the production of genetic variation, which frequently exhibits spatiotemporal Attacins (Dionne et al., 2006; Dolezal et al., 2019). The basal patterns that look as if they were shaped by selection on mode is regulated directly by dFOXO activity when induced the associated IIS traits (Fabian et al., 2012; Paaby et al., by starvation (Becker et al., 2010; Buchon et al., 2014). 2014; Kapun et al., 2016). In a recent laboratory assay, Immunity genes, including Attacins, are among the known most two foxo alleles showing opposite latitudinal clines in rapidly evolving genes and have frequently shown evidence D. melanogaster were compared on an otherwise homogeneous of local adaptation in Drosophila (Lazzaro and Clark, 2001, genetic background. The alleles showed contrasting effects 2003). on viability, size-related traits, starvation resistance, and There would be a number of mechanisms by which the fat content, whose directions were overall consistent with O7 mutation could have reduced Attacin genes’ expression. predictions from the clinal variation of the characters For example, the breakage of the invertedly transcribed AttA2 (Durmaz et al., 2019). tandem duplicates could have impaired the inducibility of one The O7 mutation could have altered Akt1 and/or foxo or the two paralogs, or their separation could have made them function via multiple non-mutually exclusive mechanisms, such lose gene expression coregulation, as might be surmised from as mutual regulatory interference, considering that they are the observations that they halted or slowed down evolving in antagonistic effectors; relocation to the sides of an immunity concert, and that AttA2b shows decreased codon bias. These gene (i.e., AttA2b) expected to be under intense purifying mechanisms could have acted synergistically with each other selection on expression (see below); and alteration of the genes’ and with those already discussed in connection with Akt1 functional neighborhood at higher-order levels of chromatin and foxo. Although this scenario could be partially offset by organization (Farré et al., 2015; McBroome et al., 2020). It the increase in basal AttA2 transcript levels that may be could be argued that the nuclear environment of the genes expected from the duplicates having been reattached to dFOXO remained basically unaltered, if the reason why they became enhancers (Becker et al., 2010), all in all, the evidence suggests involved in the rearrangement was that they already were in that (i) at its inception, O7 caused a rearrangement with close spatial proximity to each other in the nucleus. This, partial disruption of a set of functionally related loci with however, did not necessarily have to be the case, considering overlapping pleiotropic effects on immunometabolic traits. If, recent findings in yeast that rejoining of DNA break ends is in addition to these direct effects, there concurred indirect not determined by the predamage spatial proximity of the DSBs effects of linkage between locally, and given the functional (Sunder and Wilson, 2019). Be that as it may, bearing in mind relationship, likely epistatically interacting alleles warrant further that the seasonal increase of O7 occurs from early spring to investigation; and (ii) the resulting haplotype imparted a shifted Frontiers in Genetics | www.frontiersin.org 16 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance pattern of resource allocation toward reproduction at a cost to O7 may have altered fly’s immunometabolism through at immunity, compared to the OST ancestor. Such an opposing least direct effects on core immunity and metabolism genes. antagonistic pleiotropy would result in a seasonal frequency This result could help to explain the inversion’s conflicting cycle qualitatively similar to that shown by the inversions, if correlations with the seasonal and decadal climate changes, reproduction is favored from early spring to midsummer, when taking into account recent findings from microbial ecology, O7 rises (and OST wanes), and immunity from late summer to which indicate that microbial community responses to short- winter, when it wanes (and OST rises). There is ample evidence and long-term climate changes can be largely uncorrelated that the qualitative and quantitative composition of temperate (Romero-Olivares et al., 2017). Considering its large size, it bacterial communities cycles seasonally (Lazzaro et al., 2015; seems likely that O7 ’s evolution is also shaped by additional Shigyo et al., 2019). Recently, a study using D. melanogaster direct or/and indirect effects on genes other than those near from the eastern United States (Behrman et al., 2018) found its breakpoints. Further progress along this line will include a seasonal shift in immunocompetence, with the trait value development of functional tests of the identified genes on declining every spring to autumn. The shift was interpreted as inverted versus uninverted chromosome backgrounds and use resulting from relaxed selection for immune response during of the obtained assembly for building a SNP panel for O the warm season, much like what we propose here for the chromosome-wide scans of selection. We have incorporated O7 /OST inversion polymorphism. Prior data on temporal genetic the chromosome-scale sequence of O3+4+7 obtained here into variation within and between O inversions point to additional our reference genome browser9 to facilitate the further use loci that would be consistent with the seasonal cycle of O7 being of this resource. mediated by immunometabolic selection (Rodríguez-Trelles, 2003). The case of the Mpi gene encoding the key glycolytic enzyme mannose-6-phosphate isomerase (MPI) is noteworthy. From our assembly, Mpi is located 2.15 Mb outward from the DATA AVAILABILITY STATEMENT distal breakpoint of O7 , which is within the estimated region Datasets presented in this article are available at the European of the inversion-associated strong recombination–suppression Nucleotide Archive (ENA) under the project ID: PRJEB38585. effect (3.5 Mb; Pegueroles et al., 2010b). The MPI fast/slow electrophoretic polymorphism was found to be only moderately associated with the O7 /OST polymorphism. Yet (i) the magnitude of the locus-by-inversion disequilibrium cycled seasonally, and AUTHOR CONTRIBUTIONS (ii) the cycling occurred because the Fast allele increased in frequency every winter only within the O7 chromosomal class, CK, RT, and FR-T contributed to the design and implementation but not within the OST class (Rodríguez-Trelles, 2003). The of the research, to the analysis of the results, and to the writing behavior of Mpi could be in part an outcome of hitch-hiking of the manuscript. All authors contributed to the article and with other linked loci involved in seasonal adaptation. One approved the submitted version. such candidate could be the Na pumpα subunit (Atpα) gene, located only 0.13 Mb farther away from O7 than Mpi, and recently found to be under positive selection for defense against plant secondary compounds in D. subobscura (Pegueroles et al., FUNDING 2016). Still, immune elicitation in Drosophila relies upon massive This study was supported by the Spanish Ministerio de upregulation of glycolysis (Dolezal et al., 2019), which should Ciencia e Innovación grant CGL2017-89160P; and Generalitat place a strong demand on MPI activity (Shtraizent et al., 2017). de Catalunya grant 2017SGR 1379 to the Grup de Genòmica, In addition to the evidence from D. subobscura just discussed, Bioinformàtica i Biologia Evolutiva, Universitat Autònoma de Supplementary Table 2 provides additional loci found to exhibit Barcelona (Spain). CK was supported by a PIF Ph.D. fellowship seasonal variation in a genomic survey from other Drosophila, from the Universitat Autònoma de Barcelona (Spain). Note that which may be candidates for being involved in the seasonal the funding agencies were not involved in the design of the study cycling of O7 . or in any aspect of the collection, analysis and interpretation of the data or paper writing. CONCLUSION AND OUTLOOK Previous work on the spatiotemporal distribution patterns SUPPLEMENTARY MATERIAL of the inversion polymorphisms of D. subobscura indicated that O7 is driven by selective factors other than temperature The Supplementary Material for this article can be found online alone. Here, we addressed this issue using a genome-based at: https://www.frontiersin.org/articles/10.3389/fgene.2020. approach to isolate and characterize the O7 breakpoints. Our 565836/full#supplementary-material findings have general implications for current theories on the molecular mechanisms of formation of this common type of structural genomic change. Furthermore, they suggest that 9 http://dsubobscura.serveftp.com/ Frontiers in Genetics | www.frontiersin.org 17 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance REFERENCES Carvalho, C. M., Pfundt, R., King, D. A., Lindsay, S. J., Zuccherato, L. W., and Macville, M. V. (2015). Absence of heterozygosity due to template switching Aguilera, A., and Gómez-González, B. (2008). Genome instability: a mechanistic during replicative rearrangements. Am. J. Hum. Genet. 96, 555–564. doi: 10. view of its causes and consequences. Nat. Rev. Genet. 9, 204–217. doi: 10.1038/ 1016/j.ajhg.2015.01.021 nrg2268 Cer, R. Z., Bruce, K. H., Donohue, D. E., Temiz, N. A., Mudunuri, U. S., Yi, M., Alexander, J. L., Beagan, K., Orr-Weaver, T. L., and McVey, M. (2016). Multiple et al. (2012). Searching for non-B DNA-forming motifs using nBMST (non-B mechanisms contribute to double-strand break repair at rereplication forks DNA motif search tool). Curr. Protoc. Hum. Genet. 18, 11–22. doi: 10.1002/ in Drosophila follicle cells. Proc. Natl. Acad. Sci. U.S.A. 113, 13809–13814. 0471142905.hg1807s73 doi: 10.1073/pnas.1617110113 Chakraborty, M., Baldwin-Brown, J. G., Long, A. D., and Emerson, J. J. (2016). Andjelković, M., Jones, P. F., Grossniklaus, U., Cron, P., Schier, A. F., Dick, M., Contiguous and accurate de novo assembly of metazoan genomes with modest et al. (1995). Developmental regulation of expression and activity of multiple long read coverage. Nucleic Acids. Res. 44:e147. doi: 10.1093/nar/gkw654 forms of the Drosophila RAC protein kinase. J. Biol. Chem. 270, 4066–4075. Chen, H., Rangasamy, M., Tan, S. Y., Wang, H., and Siegfried, B. D. (2010). doi: 10.1074/jbc.270.8.4066 Evaluation of five methods for total DNA extraction from Western corn Anisimova, M., Gil, M., Dufayard, J. F., Dessimoz, C., and Gascuel, O. (2011). rootworm beetles. PLoS One 5:e11963. doi: 10.1371/journal.pone.0011963 Survey of branch support methods demonstrates accuracy, power, and Cheng, C., Tan, J. C., Hahn, M. W., and Besansky, N. J. (2018). Systems genetic robustness of fast likelihood-based approximation schemes. Syst. Biol. 60, analysis of inversion polymorphisms in the malaria mosquito Anopheles 685–699. doi: 10.1093/sysbio/syr041 gambiae. Proc. Natl. Acad. Sci. U.S.A. 115, E7005–E7014. doi: 10.1073/pnas. Apweiler, R., Bairoch, A., Wu, C. H., Barker, W. C., Boeckmann, B., Ferro, S., et al. 1806760115 (2004). UniProt: the universal protein knowledgebase. Nucleic Acids Res. 32, Cirera, S., Martín-Campos, J. M., Segarra, C., and Aguadé, M. (1995). Molecular D115–D119. doi: 10.1093/nar/gkw1099 characterization of the breakpoints of an inversion fixed between Drosophila Arenas, C., Zivanovic, G., and Mestres, F. (2018). Chromosomal thermal index: melanogaster and D. subobscura. Genetics 139, 321–326. a comprehensive way to integrate the thermal adaptation of Drosophila Corbett-Detig, R. B. (2016). Selection on inversion breakpoints favors proximity subobscura whole karyotype. Genome 61, 73–78. to pairing sensitive sites in Drosophila melanogaster. Genetics 204, 259–265. Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., and Cherry, J. M. doi: 10.1534/genetics.116.190389 (2000). Gene ontology: tool for the unification of biology. The Gene Ontology Corbett-Detig, R. B., and Hartl, D. L. (2012). Population genomics of inversion Consortium. Nat. Genet. 25, 25–29. doi: 10.1038/75556 polymorphisms in Drosophila melanogaster. PLoS Genet. 8:e1003056. doi: 10. Ayala, F. J., Serra, L., and Prevosti, A. (1989). A grand experiment in evolution: 1371/journal.pgen.1003056 the Drosophila subobscura colonization of the Americas. Genome 31, 246–255. Crescente, J. M., Zavallo, D., Helguera, M., and Vanzetti, L. S. (2018). doi: 10.1139/g89-042 MITE Tracker: an accurate approach to identify miniature inverted-repeat Bächli, G. (2020). TaxoDros: The Database on Taxonomy of Drosophilidae. Available transposable elements in large genomes. BMC Bioinformatics 19:348. doi: 10. online at: https://www.taxodros.uzh.ch (accessed February 21, 2020). 1186/s12859-018-2376-y Bailey, T. L., Johnson, J., Grant, C. E., and Noble, W. S. (2015). The MEME Suite. de Frutos, R. (1972). Contribution to the study of chromosomal polymorphism in Nucleic Acids Res. 43, W39–W49. doi: 10.1093/nar/gkv416 the Spanish populations of Drosophila subobscura. Genét. Ibér. 24, 123–140. Balanyà, J., Oller, J. M., Huey, R. B., Gilchrist, G. W., and Serra, L. (2006). Global de Jong, G., and Bochdanovits, Z. (2003). Latitudinal clines in Drosophila genetic change tracks global climate warming in Drosophila subobscura. Science melanogaster: body size, allozyme frequencies, inversion frequencies, and the 313, 1773–1775. doi: 10.1126/science.1131002 insulin-signalling pathway. J. Genet. 82, 207–223. doi: 10.1007/BF02715819 Bao, W., Kojima, K. K., and Kohany, O. (2015). Repbase update, a database of Delprat, A., Guillén, Y., and Ruiz, A. (2019). Computational sequence analysis of repetitive elements in eukaryotic genomes. Mob. DNA 6:11. inversion breakpoint regions in the cactophilic Drosophila mojavensis lineage. Becker, T., Loch, G., Beyer, M., Zinke, I., Aschenbrenner, A. C., Carrera, P., et al. J. Hered. 110, 102–117. doi: 10.1093/jhered/esy057 (2010). FOXO-dependent regulation of innate immune homeostasis. Nature Deng, W., Maust, B. S., Nickle, D. C., Learn, G. H., Liu, Y., Heath, L., et al. (2010). 463, 369–373. doi: 10.1038/nature08698 DIVEIN: a web server to analyze phylogenies, sequence divergence, diversity, Behrman, E. L., Howick, V. M., Kapun, M., Staubach, F., Bergland, A. O., Petrov, and informative sites. Biotechniques 48, 405–408. doi: 10.2144/000113370 D. A., et al. (2018). Rapid seasonal evolution in innate immunity of wild Dionne, M. S., Pham, L. N., Shirasu-Hiza, M., and Schneider, D. S. (2006). Akt Drosophila melanogaster. Proc. Biol. Sci. 285:20172599. doi: 10.1098/rspb.2017. and foxo dysregulation contribute to infection-induced wasting in Drosophila. 2599 Curr. Biol. 16, 1977–1985. doi: 10.1016/j.cub.2006.08.052 Benson, G. (1999). Tandem repeats finder: a program to analyze DNA sequences. Dobzhansky, T. (1947). Genetics of natural populations. XIV. A response of certain Nucleic Acids Res. 27, 573–580. doi: 10.1093/nar/27.2.573 gene arrangements in the third chromosome of Drosophila pseudoobscura to Bhandari, J., Karg, T., and Golic, K. G. (2019). Homolog-dependent repair natural selection. Genetics 32, 142–160. following dicentric chromosome breakage in Drosophila melanogaster. Genetics Dolezal, T., Krejcova, G., Bajgar, A., Nedbalova, P., and Strasser, P. (2019). 212, 615–630. doi: 10.1534/genetics.119.302247 Molecular regulations of metabolism during immune response in insects. Insect. Bracewell, R., Chatla, K., Nalley, M. J., and Bachtrog, D. (2019). Dynamic turnover Biochem. Mol. Biol. 109, 31–42. doi: 10.1016/j.ibmb.2019.04.005 of centromeres drives karyotype evolution in Drosophila. eLife 8:e49002. doi: Dolgova, O. (2013). Genetic and Phenotypic Differentiation in three Chromosomal 10.7554/eLife.49002 Arrangements of Drosophila subobscura. Ph.D. thesis, Universitat Autònoma de Brázda, V., Kolomazník, J., Lýsek, J., Hároníková, L., Coufal, J., and Št’astný, Barcelona, Barcelona. J. (2016). Palindrome analyser - A new web-based server for predicting and Durmaz, E., Rajpurohit, S., Betancourt, N., Fabian, D. K., Kapun, M., Schmidt, evaluating inverted repeats in nucleotide sequences. Biochem. Biophys. Res. P., et al. (2019). A clinal polymorphism in the insulin signaling transcription Commun. 478, 1739–1745. doi: 10.1016/j.bbrc.2016.09.015 factor foxo contributes to life-history adaptation in Drosophila. Evolution 73, Brehm, A., and Krimbas, C. B. (1988). The inversion polymorphism of Drosophila 1774–1792. doi: 10.1111/evo.13759 subobscura natural populations from Portugal. Genét. Ibér. 39, 235–248. Eilbeck, K., Lewis, S., Mungall, C., Yandell, M., Stein, L., Durbin, R., et al. (2005). Buchon, N., Silverman, N., and Cherry, S. (2014). Immunity in Drosophila The sequence ontology: a tool for the unification of genome annotations. melanogaster –from microbial recognition to whole-organism physiology. Nat. Genome Biol. 6:R44. doi: 10.1186/gb-2005-6-5-r44 Rev. Immunol. 14, 796–810. doi: 10.1038/nri3763 Ellinghaus, D., Kurtz, S., and Willhoeft, U. (2008). LTRharvest, an efficient Cáceres, M., Ranz, J. M., Barbadilla, A., Long, M., and Ruiz, A. (1999). Generation and flexible software for de novo detection of LTR retrotransposons. BMC of a widespread Drosophila inversion by a transposable element. Science 285, Bioinformatics 9:18. doi: 10.1186/1471-2105-9-18 415–418. doi: 10.1126/science.285.5426.415 Fabian, D. K., Kapun, M., Nolte, V., Kofler, R., Schmidt, P. S., Schlötterer, C., et al. Campbell, M. S., Holt, C., Moore, B., and Yandell, M. (2014). Genome annotation (2012). Genome-wide patterns of latitudinal differentiation among populations and curation using MAKER and MAKER-P. Curr. Protoc. Bioinformatics 48, of Drosophila melanogaster from North America. Mol. Ecol. 21, 4748–4769. 4.11.1–4.11.39. doi: 10.1002/0471250953.bi0411s48 doi: 10.1111/j.1365-294X.2012.05731.x Frontiers in Genetics | www.frontiersin.org 18 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance Faria, R., Johannesson, K., Butlin, R. K., and Westram, A. M. (2019). Evolving Holt, C., and Yandell, M. (2011). MAKER2: an annotation pipeline and genome- inversions. Trends Ecol. Evol. 34, 239–248. doi: 10.1016/j.tree.2018.12.005 database management tool for second-generation genome projects. BMC Farré, M., Robinson, T. J., and Ruiz-Herrera, A. (2015). An integrative breakage Bioinformatics 12:491. doi: 10.1186/1471-2105-12-491 model of genome architecture, reshuffling and evolution. Bioessays 37, 479–488. Hughes, L. (2000). Biological consequences of global warming: Is the signal already doi: 10.1002/bies.201400174 apparent? Trends Ecol. Evol. 15, 56–61. Finn, R. D., Attwood, T. K., Babbitt, P. C., Bateman, A., Bork, P., Bridge, A. J., Hughes, S. E., Miller, D. E., Miller, A. L., and Hawley, R. S. (2018). Female meiosis: et al. (2017). InterPro in 2017-beyond protein family and domain annotations. synapsis, recombination, and segregation in Drosophila melanogaster. Genetics Nucleic Acids Res. 45, D190–D199. doi: 10.1093/nar/gkw1107 208, 875–908. doi: 10.1534/genetics.117.300081 Finn, R. D., Coggill, P., Eberhardt, R. Y., Eddy, S. R., Mistry, J., Mitchell, A. L., et al. Huson, D. H., and Bryant, I. D. (2006). Application of phylogenetic networks in (2016). The Pfam protein families database: towards a more sustainable future. evolutionary studies. Mol. Biol. Evol. 23, 254–267. doi: 10.1093/molbev/msj030 Nucleic Acids Res. 44, D279–D285. doi: 10.1093/nar/gkv1344 Jones, P., Binns, D., Chang, H. Y., Fraser, M., Li, W., McAnulla, C., et al. (2014). Flatt, T., and Partridge, L. (2018). Horizons in the evolution of aging. BMC Biol. InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 16:93. doi: 10.1186/s12915-018-0562-z 1236–1240. doi: 10.1093/bioinformatics/btu031 Fontdevila, A., Zapata, C., Alvarez, G., Sanchez, L., Méndez, J., and Enriquez, I. Joyce, E. F., Paul, A., Chen, K. E., Tanneti, N., and McKim, K. S. (2012). Multiple (1983). Genetic coadaptation in the chromosomal polymorphism of Drosophila barriers to nonhomologous DNA end joining during meiosis in Drosophila. subobscura. I. Seasonal changes of gametic disequilibrium in a natural Genetics 191, 739–746. doi: 10.1534/genetics.112.140996 population. Genetics 105, 935–955. Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A., and Jermiin, Fragata, I., Lopes-Cunha, M., Bárbaro, M., Kellen, B., Lima, M., Santos, M. A., L. S. (2017). ModelFinder: fast model selection for accurate phylogenetic et al. (2014). How much can history constrain adaptive evolution? A real-time estimates. Nat. Methods 14, 587–589. doi: 10.1038/nmeth.4285 evolutionary approach of inversion polymorphisms in Drosophila subobscura. Kapun, M., Fabian, D. K., Goudet, J., and Flatt, T. (2016). Genomic evidence J. Evol. Biol. 27, 2727–2738. doi: 10.1111/jeb.12533 for adaptive inversion clines in Drosophila melanogaster. Mol. Biol. Evol. 33, Fuller, Z. L., Haynes, G. D., Richards, S., and Schaeffer, S. W. (2016). Genomics of 1317–1336. doi: 10.1093/molbev/msw016 natural populations: How differentially expressed genes shape the evolution of Kapun, M., and Flatt, T. (2018). The adaptive significance of chromosomal chromosomal inversions in Drosophila pseudoobscura. Genetics 204, 287–301. inversion polymorphisms in Drosophila melanogaster. Mol. Ecol. 28, 1263–1282. doi: 10.1534/genetics.116.191429 doi: 10.1111/mec.14871 Fuller, Z. L., Haynes, G. D., Richards, S., and Schaeffer, S. W. (2017). Genomics Karageorgiou, C., Gámez-Visairas, V., Tarrío, R., and Rodríguez-Trelles, F. (2019). of natural populations: evolutionary forces that establish and maintain gene Long-read based assembly and synteny analysis of a reference Drosophila arrangements in Drosophila pseudoobscura. Mol. Ecol. 26, 6539–6562. doi: 10. subobscura genome reveals signatures of structural evolution driven by 1111/mec.14381 inversions recombination-suppression effects. BMC Genomics 20:223. Fuller, Z. L., Koury, S. A., Phadnis, N., and Schaeffer, S. W. (2019). How Katoh, K., Rozewicki, J., and Yamada, K. D. (2019). MAFFT online service: chromosomal rearrangements shape adaptation and speciation: case studies multiple sequence alignment, interactive sequence choice and visualization. in Drosophila pseudoobscura and its sibling species Drosophila persimilis. Mol. Brief. Bioinformatics 20, 1160–1166. doi: 10.1093/bib/bbx108 Ecol. 28, 1283–1301. doi: 10.1111/mec.14923 Kaushal, S., and Freudenreich, C. H. (2019). The role of fork stalling and DNA Gienapp, P., Teplitsky, C., Alho, J. S., Mills, J. A., and Merilä, J. (2008). Climate structures in causing chromosome fragility. Genes Chromosomes Cancer 58, change and evolution: disentangling environmental and genetic responses. Mol. 270–283. doi: 10.1002/gcc.22721 Ecol. 17, 167–178. doi: 10.1111/j.1365-294X.2007.03413.x Kehrer-Sawatzki, H., Sandig, C. A., Goidts, V., and Hameister, H. (2005). Götz, W. (1965). Beitrag zur Kenntnis der Inversionen, Duplikationen und Breakpoint analysis of the pericentric inversion between chimpanzee Strukturtypen von Drosophila subobscura Coll. Z. Vererbungsl. 96, 285–296. chromosome 10 and the homologous chromosome 12 in humans. Cytogenet. doi: 10.1007/BF00896828 Genome Res. 108, 91–97. doi: 10.1159/000080806 Götz, W. (1967). Untersuchungen über den chromosomalen Kirkpatrick, M. (2010). How and why chromosome inversions evolve. PLoS Biol. Strukturpolymorphismus in kleinasiatischen und persischen Populationen 8:e1000501. doi: 10.1371/journal.pbio.1000501 von Drosophila subobscura Coll. Mol. Gen. Genet. 100, 1–38. doi: Kirkpatrick, M., and Barton, N. (2006). Chromosome inversions, local adaptation 10.1007/BF00425773 and speciation. Genetics 173, 419–434. doi: 10.1534/genetics.105.047985 Grant, C. E., Bailey, T. L., and Noble, W. S. (2011). FIMO: scanning for occurrences Kolb, J., Chuzhanova, N. A., Högel, J., Vasquez, K. M., Cooper, D. N., Bacolla, of a given motif. Bioinformatics 27, 1017–1018. doi: 10.1093/bioinformatics/ A., et al. (2009). Cruciform-forming inverted repeats appear to have mediated btr064 many of the microinversions that distinguish the human and chimpanzee Greig, D. (2007). A screen for recessive speciation genes expressed in the gametes genomes. Chromosome Res. 17, 469–483. of F1 hybrid yeast. PLoS Genet. 3:e21. doi: 10.1371/journal.pgen.0030021 Koren, S., Walenz, B. P., Berlin, K., Miller, J. R., Bergman, N. H., and Phillippy, Guindon, S., Dufayard, J. F., Lefort, V., Anisimova, M., Hordijk, W., and Gascuel, A. M. (2017). Canu: scalable and accurate long-read assembly via adaptive O. (2010). New algorithms and methods to estimate maximumlikelihood k-mer weighting and repeat separation. Genome Res. 27, 722–736. doi: 10.1101/ phylogenies: assessing the performance of PhyML 3.0. Syst. Biol. 59, 307–321. gr.215087.116 doi: 10.1093/sysbio/syq010 Koske, T., and Maynard-Smith, J. (1954). Genetics and cytology of Drosophila Hanson, M. A., and Lemaitre, B. (2020). New insights on Drosophila antimicrobial subobscura. X. The fifth linkage group. J. Genet. 52, 521–541. doi: 10.1007/ peptide function in host defense and beyond. Curr. Opin. Immunol. 62, 22–30. BF02985076 doi: 10.1016/j.coi.2019.11.008 Kramara, J., Osia, B., and Malkova, A. (2018). Break-induced replication: the Hastings, P. J., Ira, G., and Lupski, J. R. (2009). A microhomology-mediated break- where, the why, and the how. Trends Genet. 34, 518–531. doi: 10.1016/j.tig.2018. induced replication model for the origin of human copy number variation. PLoS 04.002 Genet. 5:e1000327. doi: 10.1371/journal.pgen.1000327 Krimbas, C. B. (1992). “The inversion polymorphism of Drosophila subobscura,” in Hill, T., and Betancourt, A. J. (2018). Extensive exchange of transposable elements Drosophila Inversion Polymorphism, eds C. B. Krimbas and J. R. Powell (Boca in the Drosophila pseudoobscura group. Mob. DNA 9:20. Raton, FL: CRC Press), 127–220. Hoang, D. T., Chernomor, O., von Haeseler, A., Minh, B. Q., and Vinh, L. S. (2018). Kulakovskiy, I. V., and Makeev, V. J. (2009). Discovery of DNA motifs recognized UFBoot2: improving the ultrafast bootstrap approximation. Mol. Biol. Evol. 35, by transcription factors through integration of different experimental sources. 518–522. doi: 10.1093/molbev/msx281 Biophysics 54, 667–674. doi: 10.1134/S0006350909060013 Hoffmann, A. A., and Rieseberg, L. H. (2008). Revisiting the impact of inversions Kunze-Mühl, E., and Müller, E. (1958). Weitere Untersuchungen uber in evolution: from population genetic markers to drivers of adaptive shifts and die chromosomale Struktur und die natürlichen Strukturtypen von speciation? Annu. Rev. Ecol. Evol. Syst. 39, 21–42. doi: 10.1146/annurev.ecolsys. Drosophila subobscura Coll. Chromosoma 9, 559–570. doi: 10.1007/BF0256 39.110707.173532 8093 Frontiers in Genetics | www.frontiersin.org 19 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance Kunze-Mühl, E., and Sperlich, D. (1955). Inversionen und chromosomale Orengo, D. J., Puerma, E., Cereijo, U., and Aguadé, M. (2019). The molecular Strukturtypen bei Drosophila subobscura Coll. Z. Vererbungsl. 87, 65–84. doi: genealogy of sequential overlapping inversions implies both homologous 10.1007/BF00308333 chromosomes of a heterokaryotype in an inversion origin. Sci. Rep. 9:17009. Kurtz, S., Phillippy, A., Delcher, A. L., Smoot, M., Shumway, M., Antonescu, C., Orengo, D. J., Puerma, E., Papaceit, M., Segarra, C., and Aguadé, M. A. (2015). et al. (2004). Versatile and open software for comparing large genomes. Genome Molecular perspective on a complex polymorphic inversion system with Biol. 5:R12. doi: 10.1186/gb-2004-5-2-r12 cytological evidence of multiply reused breakpoints. Heredity 114, 610–618. Lagesen, K., Hallin, P., Rødland, E. A., Staerfeldt, H. H., Rognes, T., and Ussery, doi: 10.1038/hdy.2015 D. W. (2007). RNAmmer: consistent and rapid annotation of ribosomal RNA Paaby, A. B., Bergland, A. O., Behrman, E. L., and Schmidt, P. S. (2014). A genes. Nucleic Acids Res. 35, 3100–3108. doi: 10.1093/nar/gkm160 highly pleiotropic amino acid polymorphism in the Drosophila insulin receptor Lazzaro, A., Hilfiker, D., and Zeyer, J. (2015). Structures of microbial communities contributes to life-history adaptation. Evolution 68, 3395–3409. doi: 10.1111/ in alpine soils: seasonal and elevational effects. Front. Microbiol. 6:1330. doi: evo.12546 10.3389/fmicb.2015.01330 Papaceit, M., Segarra, C., and Aguadé, M. (2012). Structure and population genetics Lazzaro, B. P., and Clark, A. G. (2001). Evidence for recurrent paralogous gene of the breakpoints of a polymorphic inversion in Drosophila subobscura. conversion and exceptional allelic divergence in the Attacin genes of Drosophila Evolution 67, 66–79. melanogaster. Genetics 159, 659–671. Parmesan, C. (2006). Ecological and evolutionary responses to recent climate Lazzaro, B. P., and Clark, A. G. (2003). Molecular population genetics of inducible change. Annu. Rev. Ecol. Evol. Syst. 37, 637–669. doi: 10.1146/annurev.ecolsys. antibacterial peptide genes in Drosophila melanogaster. Mol. Biol. Evol. 20, 37.091305.110100 914–923. doi: 10.1093/molbev/msg109 Pegueroles, C., Aquadro, C. F., Mestres, F., and Pascual, M. (2013). Gene flow and Lee, J. A., Carvalho, C. M., and Lupski, J. R. (2007). A DNA replication gene flux shape evolutionary patterns of variation in Drosophila subobscura. mechanism for generating nonrecurrent rearrangements associated with Heredity 110, 520–529. doi: 10.1038/hdy.2012.118 genomic disorders. Cell 131, 1235–1247. doi: 10.1016/j.cell.2007.11.037 Pegueroles, C., Araúz, P. A., Pascual, M., and Mestres, F. (2010a). A recombination Leigh Brown, A. J., and Ish-Horowicz, D. (1981). Evolution of the 87A and 87C survey using microsatellites: the O chromosome of Drosophila subobscura. heat-shock loci in Drosophila. Nature 290, 677–682. Genetica 138, 795–804. Lobachev, K. S., Rattray, A., and Narayanan, V. (2007). Hairpin- and cruciform- Pegueroles, C., Ferrés-Coy, A., Martí-Solano, M., Aquadro, C. F., Pascual, mediated chromosome breakage: causes and consequences in eukaryotic cells. M., and Mestres, F. (2016). Inversions and adaptation to the plant toxin Front. Biosci. 12, 4208–4220. doi: 10.2741/2381 ouabain shape DNA sequence variation within and between chromosomal Lopes, M., Foiani, M., and Sogo, J. M. (2006). Multiple mechanisms control inversions of Drosophila subobscura. Sci. Rep. 6:23754. doi: 10.1038/srep chromosome integrity after replication fork uncoupling and restart at 23754 irreparable UV lesions. Mol. Cell 21, 15–27. doi: 10.1016/j.molcel.2005.11.015 Pegueroles, C., Ordoñez, V., Mestres, F., and Pascual, M. (2010b). Recombination Loukas, M., Krimbas, C. B., and Vergini, Y. (1979). The genetics of Drosophila and selection in the maintenance of the adaptive value of inversions. J. Evol. subobscura populations. IX. Studies on linkage disequilibrium in four natural Biol. 23, 2709–2717. doi: 10.1111/j.1420-9101.2010.02136.x populations. Genetics 93, 497–523. Powell, J. R. (1997). Progress and Prospects in Evolutionary Biology: The Drosophila Lowe, T. M., and Chan, P. P. (2016). tRNAscan-SE On-line: integrating search and Model. Oxford: Oxford University Press. context for analysis of transfer RNA genes. Nucleic Acids Res. 44, W54–W57. Prazeres da Costa, O., González, J., and Ruiz, A. (2009). Cloning and sequencing doi: 10.1093/nar/gkw413 of the breakpoint regions of inversion 5g fixed in Drosophila buzzatii. Lowry, D. B., Popovic, D., Brennan, D. J., and Holeski, L. M. (2019). Mechanisms of Chromosoma 118, 349–360. doi: 10.1007/s00412-008-0201-5 a locally adaptive shift in allocation among growth, reproduction, and herbivore Prevosti, A., Ribo, G., Serra, L., Aguadé, M., Balaña, J., Monclus, M., et al. resistance in Mimulus guttatus. Evolution 73, 1168–1181. doi: 10.1111/evo. (1988). Colonization of America by Drosophila subobscura: experiment in 13699 natural populations that supports the adaptive role of chromosomal-inversion Lu, S., Wang, G., Bacolla, A., Zhao, J., Spitser, S., and Vasquez, K. M. (2015). polymorphism. Proc. Natl. Acad. Sci. U.S.A. 85, 5597–5600. doi: 10.1073/pnas. Short inverted repeats are hotspots for genetic instability: relevance to cancer 85.15.5597 genomes. Cell Rep. 10, 1674–1680. doi: 10.1016/j.celrep.2015.02.039 Puerma, E., Orengo, D. J., and Aguadé, M. (2016a). Multiple and diverse structural Mani, R. S., and Chinnaiyan, A. M. (2010). Triggers for genomic rearrangements: changes affect the breakpoint regions of polymorphic inversions across the insights into genomic, cellular and environmental influences. Nat. Rev. Genet. Drosophila genus. Sci. Rep. 6:36248. doi: 10.1038/srep36248 11, 819–829. doi: 10.1038/nrg2883 Puerma, E., Orengo, D. J., and Aguadé, M. (2016b). The origin of chromosomal Matzkin, L. M., Merritt, T. J., Zhu, C. T., and Eanes, W. F. (2005). The structure inversions as a source of segmental duplications in the Sophophora subgenus of and population genetics of the breakpoints associated with the cosmopolitan Drosophila. Sci. Rep. 6:30715. doi: 10.1038/srep30715 chromosomal inversion In(3R)Payne in Drosophila melanogaster. Genetics 170, Puerma, E., Orengo, D. J., and Aguadé, M. (2017). Inversion evolutionary rates 1143–1152. doi: 10.1534/genetics.104.038810 might limit the experimental identification of inversion breakpoints in non- Maynard-Smith, J., and Maynard-Smith, S. (1954). Genetics and cytology of model species. Sci. Rep. 7:17281. Drosophila subobscura. VIII. Heterozygosity, viability and rate of development. Puerma, E., Orengo, D. J., Cruz, F., Gómez-Garrido, J., Librado, P., Salguero, J. Genet. 52, 152–164. doi: 10.1007/BF02981496 D., et al. (2018). The high-quality genome sequence of the oceanic island McBroome, J., Liang, D., and Corbett-Detig, R. (2020). Fine-scale position effects endemic species Drosophila guanche reveals signals of adaptive evolution in shape the distribution of inversion breakpoints in Drosophila melanogaster. genes related to flight and genome stability. Genome Biol. Evol. 10, 1956–1969. Genome Biol. Evol. doi: 10.1093/gbe/evaa103 [Epub ahead of print]. doi: 10.1093/gbe/evy135 McKinney, J. A., Wang, G., Mukherjee, A., Christensen, L., Subramanian, S. H. S., Puerma, E., Orengo, D. J., Salguero, D., Papaceit, M., Segarra, C., and Aguadé, M. Zhao, J., et al. (2020). Distinct DNA repair pathways cause genomic instability (2014). Characterization of the breakpoints of a polymorphic inversion complex at alternative DNA structures. Nat. Commun. 11:236. detects strict and broad breakpoint reuse at the molecular level. Mol. Biol. Evol. Menozzi, P., and Krimbas, C. B. (1992). The inversion polymorphism of 31, 2331–2341. doi: 10.1093/molbev/msu177 D. subobscura revisited: synthetic maps of gene arrangement frequencies and Puig-Giribets, M., Guerreiro, M. P. G., Santos, M., Ayala, F. J., Tarrío, R., their interpretation. J. Evol. Biol. 5, 625–641. doi: 10.1046/j.1420-9101.1992. and Rodríguez-Trelles, F. (2019). Chromosomal inversions promote genomic 5040625.x islands of concerted evolution of Hsp70 genes in the Drosophila subobscura Messer, P. W., Ellner, S. P., and Hairston, N. G. Jr. (2016). Can population genetics species subgroup. Mol. Ecol. 28, 1316–1332. doi: 10.1111/mec.14511 adapt to rapid evolution? Trends Genet. 32, 408–418. doi: 10.1016/j.tig.2016.04. Ranz, J. M., Maurin, D., Chan, Y. S., von Grotthuss, M., Hillier, L. W., Roote, J., et al. 005 (2007). Principles of genome evolution in the Drosophila melanogaster species Nguyen, L.-T., Schmidt, H. A., von Haeseler, A., and Minh, B. Q. (2015). IQ-TREE: group. PLoS Biol. 5:e152. doi: 10.1371/journal.pbio.0050152 a fast and effective stochastic algorithm for estimating maximum likelihood Reese, M. G. (2001). Application of a time-delay neural network to promoter phylogenies. Mol. Biol. Evol. 32, 268–274. doi: 10.1093/molbev/msu300 annotation in the Drosophila melanogaster genome. Comput. Chem. 26, 51–56. Frontiers in Genetics | www.frontiersin.org 20 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance Regan, J. C., Froy, H., Walling, C. A., Moatt, J. P., and Nussey, D. H. (2020). Dietary Seppey, M., Manni, M., and Zdobnov, E. M. (2019). BUSCO: assessing genome restriction and insulin-like signalling pathways as adaptive plasticity: a synthesis assembly and annotation completeness. Methods Mol. Biol. 1962, 227–245. and re-evaluation. Funct. Ecol. 34, 107–128. doi: 10.1111/1365-2435.13418 doi: 10.1007/978-1-4939-9173-0_14 Rego, C., Balanyà, J., Fragata, I., Matos, M., Rezende, E. L., and Santos, M. Shigyo, N., Umeki, K., and Hirao, T. (2019). Seasonal dynamics of soil fungal (2010). Clinal patterns of chromosomal inversion polymorphisms in Drosophila and bacterial communities in cool-temperate montane forests. Front. Microbiol. subobscura are partly associated with thermal preferences and heat stress 10:1944. doi: 10.3389/fmicb.2019.01944 resistance. Evolution 64, 385–397. doi: 10.1111/j.1558-5646.2009.00835.x Shtraizent, N., DeRossi, C., Nayar, S., Sachidanandam, R., Katz, L. S., Rezende, E. L., Balanyà, J., Rodríguez-Trelles, F., Rego, C., Fragata, I., Matos, Prince, A., et al. (2017). MPI depletion enhances O-GlcNAcylation of M., et al. (2010). Climate change and chromosomal inversions in Drosophila p53 and suppresses the Warburg effect. eLife 6:e22477. doi: 10.7554/eLife. subobscura. Clim. Res. 43, 103–114. doi: 10.1007/s10709-018-0035-x 22477 Richards, S., Liu, Y., Bettencourt, B. R., Hradecky, P., Letovsky, S., and Nielsen, Slade, J. D., and Staveley, B. E. (2016). Enhanced survival of Drosophila Akt1 R. (2005). Comparative genome sequencing of Drosophila pseudoobscura: hypomorphs during amino-acid starvation requires foxo. Genome 59, 87–93. chromosomal, gene, and cis-element evolution. Genome Res. 15, 1–18. doi: Smit, A. F. A., Hubley, R., and Green, P. (2013/2015). RepeatMasker Open- 10.1101/gr.3059305 4.0. Available online at: http://www.repeatmasker.org (accessed November 15, Rivera, J., Keränen, S. V. E., Gallo, S. M., and Halfon, M. S. (2019). REDfly: the 2019). transcriptional regulatory element database for Drosophila. Nucleic Acids Res. Smith, C., Llorente, B., and Symington, L. (2007). Template switching during 47, D828–D834. doi: 10.1093/nar/gky957 break-induced replication. Nature 447, 102–105. doi: 10.1038/nature05723 Rodríguez-Trelles, F. (2003). Seasonal cycles of allozyme-by-chromosomal- Soderlund, C., Bomhoff, M., and Nelson, W. M. (2011). SyMAP v3.4: a turnkey inversion gametic disequilibrium in Drosophila subobscura. Evolution 57, 839– synteny system with application to plant genomes. Nucleic Acids Res. 39:e68. 848. doi: 10.1111/j.0014-3820.2003.tb00295.x doi: 10.1093/nar/gkr123 Rodríguez-Trelles, F., Alvarez, G., and Zapata, C. (1996). Time-series analysis of Sogo, J. M., Lopes, M., and Foiani, M. (2002). Fork reversal and ssDNA seasonal changes of the O inversion polymorphism of Drosophila subobscura. accumulation at stalled replication forks owing to checkpoint defects. Science Genetics 142, 179–187. 297, 599–602. doi: 10.1126/science.1074023 Rodríguez-Trelles, F., and Rodríguez, M. A. (1998). Rapid micro-evolution and loss Solé, E., Balanyà, J., Sperlich, D., and Serra, L. (2002). Long-term changes of chromosomal diversity in Drosophila in response to climate warming. Evol. in the chromosomal inversion polymorphism of Drosophila subobscura. I. Ecol. 12, 829–838. doi: 10.1023/A:1006546616462 Mediterranean populations from southwestern Europe. Evolution 56, 830–835. Rodríguez-Trelles, F., and Rodríguez, M. A. (2007). Comment on ‘Global genetic doi: 10.1111/j.0014-3820.2002.tb01393.x change tracks global climate warming in Drosophila subobscura’. Science Sperlich, D., Feuerbach-Mravlag, H., Lange, P., Michaelidis, A., and Pentzos- 315:1497. doi: 10.1126/science.1136298 Daponte, A. (1977). Genetic load and viability distribution in central and Rodríguez-Trelles, F., and Rodríguez, M. A. (2010). Measuring evolutionary marginal populations of Drosophila subobscura. Genetics 86, 835–848. responses to global warming: cautionary lessons from Drosophila. Insect Sun, X., Qun, Y., and Xia, X. (2013). An improved implementation of effective Conserv. Divers. 3, 44–50. doi: 10.1111/j.1752-4598.2009.00071.x number of codons (Nc). Mol. Biol. Evol. 30, 191–196. doi: 10.1093/molbev/ Rodríguez-Trelles, F., Tarrío, R., and Santos, M. (2013). Genome-wide evolutionary mss201 response to a heat wave in Drosophila. Biol. Lett. 9:20130228. doi: 10.1098/rsbl. Sunder, S., and Wilson, E. T. (2019). Frequency of DNA end joining in trans is 2013.0228 not determined by the predamage spatial proximity of double-strand breaks Romero-Olivares, A. L., Allison, S. D., and Treseder, K. K. (2017). Soil microbes in yeast. Proc. Natl. Acad. Sci. U.S.A. 116, 9481–9490. doi: 10.1073/pnas. and their response to experimental warming over time: a meta-analysis of field 1818595116 studies. Soil Biol. Biochem. 107, 32–40. doi: 10.1016/j.soilbio.2016.12.026 Tajima, F. (1993). Simple methods for testing the molecular evolutionary clock Said, I., Byrne, A., Serrano, V., Cardeno, C., Vollmers, C., and Corbett-Detig, R. B. hypothesis. Genetics 135, 599–607. (2018). Linked genetic variation and not genome structure causes widespread Teleman, A. A. (2010). Molecular mechanisms of metabolic regulation by insulin differential expression associated with chromosomal inversions. Proc. Natl. in Drosophila. Biochem. J. 425, 13–26. Acad. Sci. U.S.A. 115, 5492–5497. doi: 10.1073/pnas.1721275115 The Gene Ontology Consortium. (2017). Expansion of the gene ontology Sakofsky, C. J., Ayyar, S., Deem, A. K., Chung, W. H., Ira, G., and Malkova, knowledgebase and resources. Nucleic Acids Res. 45, D331–D338. doi: 10.1093/ A. (2015). Translesion polymerases drive microhomology-mediated break- nar/gkw1108 induced replication leading to complex chromosomal rearrangements. Mol. Cell Tremblay-Belzile, S., Lepage, É., Zampini, É., and Brisson, N. (2015). Short- 60, 860–872. doi: 10.1016/j.molcel.2015.10.041 range inversions: rethinking organelle genome stability: template switching Santos, J., Serra, L., Solé, E., and Pascual, M. (2010). FISH mapping of microsatellite events during DNA replication destabilize organelle genomes. Bioessays 37, loci from Drosophila subobscura and its comparison to related species. 1086–1094. doi: 10.1002/bies.201500064 Chromosome Res. 18, 213–226. Voineagu, I., Narayanan, V., Lobachev, K. S., and Mirkin, S. M. (2008). Replication Santos, M., Céspedes, W., Balanyà, J., Trotta, V., Calboli, F. C., Fontdevila, A., stalling at unstable inverted repeats: interplay between DNA hairpins and fork et al. (2005). Temperature-related genetic changes in laboratory populations of stabilizing proteins. Proc. Natl. Acad. Sci. U.S.A. 105, 9936–9941. doi: 10.1073/ Drosophila subobscura: evidence against simple climatic-based explanations for pnas.0804510105 latitudinal clines. Am. Nat. 165, 258–273. doi: 10.1086/427093 Walker, B. J., Abeel, T., Shea, T., Priest, M., Abouelliel, A., Sakthikumar, S., et al. Savage, J. R. (1976). Classification and relationships of induced chromosomal (2014). Pilon: an integrated tool for comprehensive microbial variant detection structual changes. J. Med. Genet. 13, 103–122. doi: 10.1136/jmg. and genome assembly improvement. PLoS One 9:e112963. doi: 10.1371/journal. 13.2.103 pone.0112963 Schwenke, R. A., Lazzaro, B. P., and Wolfner, M. F. (2016). Reproduction- Wang, G., Christensen, L. A., and Vasquez, K. M. (2006). Z-DNA-forming immunity trade-offs in insects. Annu. Rev. Entomol. 61, 239–256. sequences generate large-scale deletions in mammalian cells. Proc. Natl. Acad. Scully, R., Panday, A., Elango, R., and Willis, N. A. (2019). DNA double-strand Sci. U.S.A. 103, 2677–2682. doi: 10.1073/pnas.0511084103 break repair-pathway choice in somatic mammalian cells. Nat. Rev. Mol. Cell Wang, G., and Vasquez, K. M. (2006). Non-B DNA structure-induced genetic Biol. 20, 698–714. instability. Mutat. Res. 598, 103–119. Sebastián, A., and Contreras-Moreira, B. (2014). FootprintDB: a database of Wasserman, M. (1968). Recombination-induced chromosomal heterosis. Genetics transcription factors with annotated cis elements and binding interfaces. 58, 125–139. Bioinformatics 30, 258–265. doi: 10.1093/bioinformatics/btt663 Weirauch, M. T., Yang, A., Albu, M., Cote, A. G., Montenegro-Montero, A., Senger, K., Armstrong, G. W., Rowell, W. J., Kwan, J. M., Markstein, M., and Drewe, P., et al. (2014). Determination and inference of eukaryotic transcription Levine, M. (2004). Immunity regulatory DNAs share common organizational factor sequence specificity. Cell 158, 1431–1443. doi: 10.1016/j.cell.2014. features in Drosophila. Mol. Cell 13, 19–32. 08.009 Frontiers in Genetics | www.frontiersin.org 21 October 2020 | Volume 11 | Article 565836 Karageorgiou et al. Inversion O7 Origin and Maintenance Wellenreuther, M., and Bernatchez, L. (2018). Eco-evolutionary genomics of Zouros, E., Krimbas, C. B., Tsakas, S., and Loukas, M. (1974). Genic versus chromosomal inversions. Trends Ecol. Evol. 33, 427–440. doi: 10.1016/j.tree. chromosomal variation in natural populations of D. subobscura. Genetics 78, 2018.04.002 1223–1244. Wesley, C. S., and Eanes, W. F. (1994). Isolation and analysis of the breakpoint sequences of chromosome inversion In(3L)Payne in Drosophila melanogaster. Conflict of Interest: The authors declare that the research was conducted in the Proc. Natl. Acad. Sci. U.S.A. 91, 3132–3136. doi: 10.1073/pnas.91.8.3132 absence of any commercial or financial relationships that could be construed as a Zhang, F., Khajavi, M., Connolly, A. M., Towne, C. F., Batish, S. D., and Lupski, potential conflict of interest. J. R. (2009). The DNA replication FoSTeS/MMBIR mechanism can generate genomic, genic and exonic complex rearrangements in humans. Nat. Genet. 41, Copyright © 2020 Karageorgiou, Tarrío and Rodríguez-Trelles. This is an open-access 849–853. doi: 10.1038/ng.399 article distributed under the terms of the Creative Commons Attribution License Zhao, J., Bacolla, A., Wang, G., and Vasquez, K. M. (2010). Non-B DNA structure- (CC BY). The use, distribution or reproduction in other forums is permitted, provided induced genetic instability and evolution. Cell. Mol. Life Sci. 67, 43–62. the original author(s) and the copyright owner(s) are credited and that the original Zollinger, E. (1950). Ein strukturell homozygoter Stamm von Drosophila publication in this journal is cited, in accordance with accepted academic practice. No subobscura aus einer Wild-population. Arch. Klaus-Stiftg 25, 33–35. use, distribution or reproduction is permitted which does not comply with these terms. Frontiers in Genetics | www.frontiersin.org 22 October 2020 | Volume 11 | Article 565836