Differences in orthodenticle expression promote ommatidial size 1 variation between Drosophila species 2

28 The compound eyes of insects exhibit extensive variation in ommatidia number and size, 29 which affects how they see and underlies adaptations in their vision to different 30 environments and lifestyles. However, very little is known about the genetic and 31 developmental bases underlying differences in compound eye size. We previously showed 32 that the larger eyes of Drosophila mauritiana compared to D. simulans is caused by 33 differences in ommatidia size rather than number. Furthermore, we identified an X-linked 34 chromosomal region in D. mauritiana that results in larger eyes when introgressed into D. 35 simulans . Here, we used a combination of fine-scale mapping and gene expression 36 analysis to further investigate positional candidate genes on the X chromosome. We found 37 that orthodenticle is expressed earlier in D. mauritiana than in D. simulans during 38 ommatidial maturation in third instar larvae, and we further show that this gene is required 39 for the correct organisation and size of ommatidia in D. melanogaster . Using ATAC-seq, 40 we have identified several candidate eye enhancers of otd as well as potential direct 41 targets of this transcription factor that are differentially expressed between D. mauritiana 42 and D. simulans . Taken together, our results suggest that differential timing of otd 43 expression contributes to natural variation in ommatidia size between D. mauritiana and D. 44 simulans, which provides new insights into the mechanisms underlying the regulation and 45 evolution of compound eye size in insects. 46


Introduction
Understanding the genetic and genomic basis of phenotypic diversity is one of the central themes of evolutionary biology.While the causative genes and even mutations have been identified underlying evolutionary changes in a growing list of phenotypes (e. g. (Arif et al., 2013;Arnoult et al., 2013;Hagen et al., 2019;Klaassen et al., 2018;Ramaekers et al., 2019;Santos et al., 2017) and see (Courtier-Orgogozo et al., 2019) for a more comprehensive list) we still know relatively little about the genetic basis for the evolution of organ size.Identifying such genes will not only broaden our understanding of morphological change but provide further insights into the mechanisms underlying the control of organ size.
Insects exhibit remarkable variation in the size and shape of their compound eyes, which has allowed these animals to adapt to different environments and lifestyles (Land and Nilsson, 2012).This variation greatly affects optical parameters and visual sensation, such as the detection of different intensities, polarization and wavelengths of light to varying degrees of contrast sensitivity and acuity (Land and Nilsson, 2012).Compound eyes vary in the size and/or number of ommatidia: wider ommatidia capture more light, which can increase contrast sensitivity; however, larger interommatidial angles can lead to decreased acuity (Land, 1997;Land and Nilsson, 2012).Conversely, having many small ommatidia with narrow interommatidial angles can enhance acuity, but this may decrease contrast sensitivity (Currea et al., 2018;Palavalli-Nettimi and Theobald, 2020;Warrant, 1999Warrant, , 2006)).Differences in ommatidia number and size, as well as trade-offs between these structural features of compound eyes, have been described for a range of different insects (Duncan et al., 2020;Gonzalez-Bellido et al., 2011;Horridge, 1977;Posnien et al., 2012;Wakakuwa et al., 2007).Furthermore, variation in ommatidia size across the eye within species is also widely documented (Land, 1989).This size variation demonstrates areas of regional specialisation, where different visual tasks are performed by different parts of the eye.For example, killer flies, Coenosia sp., have evolved wider, flattened, anterior ommatidia to maximise contrast sensitivity and acuity as an adaptation to hunting (Gonzalez-Bellido et al., 2011).A number of studies have also found extensive variation in eye size within and between closely related species of Drosophila, caused by differences in ommatidia number and/or ommatidia area (Arif et al., 2013;Buchberger et al., 2021;Gaspar et al., 2020;Hilbrant et al., 2014;Keesey et al., 2019;Norry and Gomez, 2017;Posnien et al., 2012;Ramaekers et al., 2019;Reis et al., 2020).Despite the pervasive variation in eye morphology and detailed knowledge about eye development in Drosophila .CC-BY 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is The copyright holder for this preprint this version posted March 17, 2021.; https://doi.org/10.1101/2021.03.17.435774 doi: bioRxiv preprint (Casares and Almudi, 2016;Casares and McGregor, 2021;Kumar, 2018), little is known about the genetic and developmental bases for variation in eye size even among Drosophila species with a few exceptions (e.g.Ramaekers et al., 2019).
We previously showed that D. mauritiana has larger eyes than D. simulans due to larger ommatidia rather than an increase in ommatidia number (Arif et al., 2013;Posnien et al., 2012).Quantitative trait loci (QTL) mapping of this difference identified a large-effect QTL that explains 33% of the species difference.Introgression of this X-linked region from D. mauritiana into D. simulans increased the eye size and ommatidial size of the latter species (Arif et al., 2013).Here, we combine higher resolution mapping of this previously characterised X-linked QTL, with transcriptomic analysis of eye imaginal discs of D. simulans and D. mauritiana, to identify positional candidate genes that are differentially expressed in the developing ommatidia between these two species.We then carry out ATAC-seq to compare putative regulatory regions of the candidate gene orthodenticle (otd) that may contribute to differences in ommatidia diameter between D. mauritiana and D. simulans.Our results suggest that differential regulation of otd results in earlier expression of this homeobox gene in D. mauritiana compared to D. simulans.We hypothesise that this heterochrony in otd expression and consequently longer exposure to this transcription factor (TF) in maturing ommatidia in D. mauritiana contributes to the development of larger ommatidia in this species.

Enlarged ommatidia in D. mauritiana
We previously found that the larger eyes of D. mauritiana compared to those of D. simulans are caused by wider diameter of central ommatidia in the former species (Arif et al., 2013;Posnien et al., 2012).To examine whether this phenotypic difference is prevalent in all ommatidia across the eye, we imaged the eyes of these two species using synchrotron radiation micro CT (SRµCT) and measured the facet diameter of ommatidia in different regions of the eye using a 3D reconstruction of each species (Fig. 1).We corroborated that while the number of ommatidia is similar between D. mauritiana and D.
simulans the former has larger facets.This is consistent across anterior, central and posterior facets but is particularly pronounced in the antero-ventral region of the eye (Fig. 1 and Suppl.Table 1) (Arif et al., 2013;Posnien et al., 2012)

Differentially expressed genes in a candidate region on the X chromosome
Previously we detected a QTL region located between 2.6 Mb and 8 Mb on the X chromosome, which is responsible for 33% of the difference in ommatidia size (Arif et al., 2013).Furthermore, introgression of approximately 8.3 Mb of this X-linked region (between the yellow (y) and the vermillion (v) loci) from D. mauritiana TAM16 into D. simulans YVF significantly increased the eye size of the latter, consistent with the direction of the species difference (Arif et al., 2013).Further analysis of recombinant males with breakpoints within the introgressed region revealed significant genotype-phenotype associations towards the distal end of the introgressed region near marker v, providing a conservative interval of about 2 Mb wherein the X QTL is likely to reside (Fig. 2a).To map the candidate region to higher resolution we generated introgression lines with breakpoints in the 2 Mb interval and compared eye area and central ommatidia diameter of yf male progeny (with some D. mauritiana DNA in the 2Mb interval) to that of their yvf sibling males (i.e., without D. mauritiana DNA).We found that yf males had significantly larger eye size than their yvf siblings in introgression lines IL9.1a (one tailed t=1.80, df=11, p=.026) IL9.1b (one tailed t=1.80, df=11, p<.001) and IL9.2 (one tailed t=1.80, df=11, p=.035) but ommatidia diameter was only significantly different for IL9.1 (one tailed t=1.80, df=11, p=.014) and b (one tailed t=1.80, df=11, p=.005) (Fig. 2a).Ommatidia number and body size did not differ between yf males and their respective yvf sibling males for any of the IL lines (Suppl. Table 2).These data suggest the candidate QTL is located in a maximum region of just over 662 kb (ChrX: 7,725,387,618 in D. simulans)   This mapped region contains 62 protein coding positional candidate genes.To assay which positional candidates are actually expressed during the generation of .CC-BY 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is The copyright holder for this preprint this version posted March 17, 2021.; https://doi.org/10.1101/2021.03.17.435774 doi: bioRxiv preprint ommatidia, we performed RNA-seq experiments on the eye-antennal imaginal discs of 3rd instar larvae (L3).We extracted RNA from D. mauritiana and D. simulans eye-antennal discs at three different developmental points: at 72 hours after egg laying (AEL; late L2, at the onset of the morphogenetic furrow (MF)), when the eye primordium is proliferating and specification of the ommatidial cells has not yet started; at 96 h AEL stage (mid L3) when the MF has moved half way through the eye disc) and the most posterior ommatidia are already determined, and finally, at 120 h AEL (late L3), when the larvae are about to pupariate and most ommatidia are already determined whilst their final size, structure and shape are being arranged (Torres-Oliva et al., 2016;Torres-Oliva et al., 2018).
Comparison of the RNA-seq data among these three developmental timepoints showed that transcriptomes of 72 h AEL eye imaginal discs were the most different in comparison to transcriptomes from both 96 h AEL and 120 h AEL for both species (Suppl.Fig. 1).This reflects the distinctive processes that are occurring at these developmental stages (Torres-Oliva et al., 2018).We next focused on the expression of genes located within the mapped 0.66 Mb X-linked region and, in particular, on expression differences at 120 h AEL, when at least the most posterior ommatidia are acquiring their final size.Of the 62 genes located in this region, 49 were expressed in at least one of our RNA-seq datasets and only eight of these genes were differentially expressed between the eye discs of these two species at this timepoint (Suppl.Table 3): spirit, otd and Ppt1 showed higher expression in D. mauritiana, whereas CG1632, Es2, Sptr, CG12112 and CG1885 were more highly expressed in D. simulans (Fig. 2b).
We next performed in situ hybridization experiments of these eight candidate genes to investigate if they are expressed in the eye field where the ommatidia are being assembled.These assays were carried out in both D. mauritiana and D. simulans, which allowed us to determine whether the differences in expression levels observed in the RNAseq datasets are related to differences in spatial expression (Fig. 2c).Sptr, CG12112 and spirit had no detectable expression in the relevant region posterior to the MF (Fig. 2c).
Ppt1 and CG1885 were expressed both anterior to and immediately posterior to the MF.CG1632 and Es2 were ubiquitously expressed in the eye disc, with no clear regional differences.Finally, otd was expressed in the ocellar region of the disc and in the most posterior region of the eye field, where the ommatidia are already determined and are being assembled (Vandendries et al., 1996).Otd is already expressed in several rows of most posterior ommatidia of D. mauritiana eye discs, whereas, otd expression is undetectable in the most posterior regions of the eye discs of D. simulans (Fig. 2c).These results were consistent with our differential expression (DE) analysis, as there appeared to .CC-BY 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is The copyright holder for this preprint this version posted March 17, 2021.; https://doi.org/10.1101/2021.03.17.435774 doi: bioRxiv preprint be detected qualitative differences in expression levels for most of the genes investigated.
Taken together, these results showed that otd is the only differentially expressed positional candidate gene that is expressed in maturing ommatidia (Fig. 2c).

Differences in otd gene expression during eye development between D. simulans and D. mauritiana
Our results suggested that otd transcription in the maturing ommatidia initiates earlier in D. mauritiana than in D. simulans eye discs (Fig. 2c).To investigate this further, we performed additional in situ hybridizations at 110 hAEL to compare the onset of otd expression in the developing ommatidia of these two species.At this developmental stage, we found that otd is already expressed in D. mauritiana eye discs, whereas there was no detectable expression in D. simulans discs (Fig. 3a-d).To confirm this heterochrony in otd expression we performed immunostainings against Otd protein in developing eye discs (Fig. 3e, f).We measured the number of ommatidial rows that were already specified (i.e. with positive Elav staining) as a proxy of developmental stage and then which of these ommatidial rows showed Otd expression.We observed that D. mauritiana discs displayed more Otd-positive ommatidia than D. simulans discs at the same stage (Elav positive ommatidia rows, Fig. 3e-g, Suppl.Table 4, F 1,47 = 30.3, p-value=1.48x 10 -06 ).Thus, D. mauritiana cells in maturing ommatidia are exposed to the action of the TF Otd for longer since the expression of this protein extends into the pupal stage of both species (Suppl.otd is required for the correct arrangement and size of ommatidia in Drosophila It was previously shown in D. melanogaster that otd is expressed in the photoreceptor cells and is required during early pupal stages for morphogenesis of the rhabdomeres, and subsequently rhodopsin expression, as well as synaptic-column and layer targeting of these cells (Fichelson et al., 2012;Mencarelli and Pichaud, 2015;Vandendries et al., 1996).We carried out further analysis of otd function during eye development using RNAi knockdown and by generating mitotic clones of homozygous otd mutant cells.Decreasing otd mRNA by overexpressing otd-miRNA construct (Wang et al., 2010), resulted in defects in the final ommatidia organisation that were rescued by adding a copy of UAS-otd (Suppl.Fig. 3a-c).Loss of otd in clones resulted in disorganised ommatidia with perturbed shapes and sizes -often smaller than those of ommatidia of controls (Suppl.Fig. 3d).These results show that otd expression in the photoreceptor cells of maturing ommatidia is required for the proper regulation of ommatidial organisation and size.To investigate the regulation of otd in the developing eyes of D. simulans and D. mauritiana further, we performed ATAC-seq (Buchberger et al., 2021;Buenrostro et al., 2013;Kittelmann et al., 2018) on D. simulans and D. mauritiana eye imaginal discs at 96 and 120 h AEL.We mapped our datasets against both genomes in order to detect common, differentially accessible and species-specific regulatory regions (Fig. 4a).The ATAC-seq peak calling of the four datasets (two developmental stages and two species) revealed a total of 22 peaks in the otd locus, all of which were located within alignable orthologous regions in the two species (Fig. 4a).Four of these peaks showed significant differences in accessibility between D. mauritiana and D. simulans: peak 6 (D. sim chrX: 8,100,587-8,100,808, padj = 0.00155) and peak 11 (D.sim chrX: 8,107,881-8,108,402, padj = 0.00976) in the 3rd and 1st introns of otd, respectively, and peaks 17 (D. sim chrX: 8,125,765 -8,126,410, padj = 0.0418) and 19 (D. sim chrX: 8,129,876 -8,131,170, padj = 0.0317) located upstream of otd (Fig. 4).

Differences in chromatin accessibility
We aligned the orthologous sequences of these differential peaks and found that each of them contained several sequence variants (peak 6: 8 SNPs and 2 small indels; peak 11: 8 SNPs and 5 small indels; peak 17: 4 SNPs and 1 small indel and peak 19: 29 SNPs and 4 small indels, Suppl.Fig. 5).Furthermore, we found a fifth peak (peak 16) upstream of otd, which did not show differential accessibility within the alignable regions, but high accessibility specifically in two D. mauritiana stretches that are disrupted by two insertions of 55 and 106 nucleotides respectively in D. simulans (Fig. 4b-d).
We next looked for transcription factor binding motifs (TFBMs) corresponding to the variant sequences within the four differential accessible peaks and in the two D. mauritiana-specific regions in peak 16.We predicted 164 putative TFBM in D. simulans and D. mauritiana in Peak 6 using JASPAR with a threshold of 85%.This included predicted TFBMs for Cut (Ct), Odd-paired (Opa), Optix, Sine oculis (So) and the Iroquois complex TFs (Araucan (Ara), Caupolican (Caup) and Mirror), which all have detectable expression in our RNA-seq datasets (Suppl.Table 3) and are known to be involved in eye-  5) and 60 were D. simulans-specific (Trl, Dll, Optix, Iroquois Complex, Dref, BEAF-32, Ems, Cut, Deaf1, E5, Ovo, Sd, Bgb::run).Overall, we identified several regions in the otd locus that may function as new eye enhancers for this gene in D. mauritiana and D. simulans, and that their differential activity could be responsible for the temporal differences in the onset of otd activation between these two species.

Differences in otd targets during eye development between D. simulans and D. mauritiana
Next, we investigated whether differences in the onset of expression of otd between D. simulans and D. mauritiana promoted further changes in its Gene Regulatory Networks (GRNs) that ultimately may be responsible for the differences in the final size of ommatidia.To this end, we called open chromatin peaks in each sample and searched for the Otd-binding motif in these accessible regions of genes expressed during eye development.Based on this analysis, we found 1,148 putative Otd target genes.We next examined which of these accessible chromatin peaks were associated with genes that were differentially expressed in our transcriptome datasets.We found that peaks associated with 161 of the 1330 genes that are upregulated in D. mauritiana contained Otd binding motifs, and 111 out of 1249 genes upregulated in D. simulans had associated peaks containing Otd binding motifs (Fig. 5a, c, Suppl.Table 6).We then performed Gene Ontology (GO) term enrichment analysis for these differentially expressed genes with accessible chromatin containing Otd binding motifs.
The D. mauritiana dataset exhibited enrichment in terms related to developmental processes, cell morphogenesis, axonogenesis, regulation of differentiation or growth, among others (Fig. 5b).By contrast, genes that were upregulated in D. simulans with associated Otd-peaks were enriched in terms such as the MAP kinase network, signalling by interleukins and cellular response to insulin stimulus (Fig. 5d).

Ommatidia size differences related to Rh3 expression
Given that Otd directly regulates the expression of rh3 (Tahayato et al., 2003) and we have previously shown higher levels of rh3 expression in adult eyes in D. mauritiana (Hilbrant et al., 2014;Posnien et al., 2012), we tested if the expression of Rh3 in ommatidia may have an impact on the ommatidia diameter.We used immunohistology to detect Rh3 + ommatidia and applied confocal microscopy as well as 3D reconstructions to measure ommatidia in a 10x10 quadrant in the central eye region (Fig. 6a, b).This analysis showed that Rh3 + ommatidia were wider in both species, while this difference was more pronounced in D. mauritiana (Fig. 6d).Note that we did not find differences in the ratio of Rh3 + /Rh3 -ommatidia between species (Fig. 6c).This data suggested exposing maturing ommatidia longer to the action of Otd may quantitatively influence the Rh3 content, rather than the ratio of ommatidia subtypes.7.

Discussion
In Drosophila, ommatidia are produced in a posterior-to-anterior hexagonal pattern in the wake of the morphogenetic furrow (MF), which passes the presumptive retinal field during L3 to trigger photoreceptor cell specification and differentiation, and subsequently the formation of cone cells and other ommatidial cells (reviewed in (Casares and Almudi, 2016;Casares and McGregor, 2021;Gaspar et al., 2019;Kumar, 2018)).Then in pupal stages, the cone cells are placed above the photoreceptor cells and this, in combination with apoptosis and further rearrangement of pigment and other cells, produces the final number and arrangement of ommatidial cells.Final ommatidia size is specified by about 40 hours after puparium formation (hAPF), the lens is secreted from 60 hAPF, and the rhodopsins are expressed from 96 hAPF (Cagan and Ready, 1989;Earl and Britt, 2006;Kim et al., 2016).
While much is known about the specification and differentiation of ommatidia, very little is known about the regulation and evolution of their size.To investigate the genetic basis of the difference in ommatidia size between D. mauritiana and D. simulans, we carried out high-resolution introgression mapping of a previously identified X linked QTL previously mapped that explains about 33% of the difference in eye size between these two species (Arif et al., 2013).We identified eight positional candidate genes in this region that differed in expression between the developing eye-antennal discs of D. mauritiana and D. simulans.Our analysis of the spatial expression of these eight genes strongly suggests otd as the best candidate gene in this region for the underlying difference in ommatidia size and thus overall eye size between these species.
otd/Otx genes play several important roles during eye development in both invertebrates and vertebrates (Ragge et al., 2005;Ranade et al., 2008;Sen et al., 2013;Tahayato et al., 2003;Vandendries et al., 1996).During eye development, Otd regulates genes for cell adhesion and cytoskeletal organisation and this is essential for the correct development of the photoreceptor cells and ommatidia maturation (Fichelson et al., 2012;Ranade et al., 2008).Mutations in otd perturb morphogenesis of the photoreceptor cells which affects the spacing of the more apical cone cells (Fichelson et al., 2012;Vandendries et al., 1996).Intriguingly, the removal of photoreceptor cells changes ommatidia size (Miller and Cagan, 1998).We propose that although otd is not expressed in the lens-secreting cone cells, it indirectly affects the organisation of these cells and thus ommatidia size through regulating the maturation and organisation of the underlying photoreceptor cells.We have shown that knockdown or loss of otd in D. melanogaster perturbs ommatidia size specification, but it remains to be directly tested if variation in the expression of this gene underlies larger and smaller ommatidia in D. mauritiana and D.
simulans respectively and if otd contributes to the observed variation in ommatidia size in different regions of the eye.
Changes in developmental timing, or heterochrony, have played an essential role in the evolution of morphologies in multiple taxa (Alberch and Alberch, 1981;Alberch et al., 1979;Gould, 1977;McKinney, 1988).Classically, the term heterochrony has been used to refer to differences in the timing of developmental events and several examples of heterochrony have been described (Briscoe and Small, 2015;Ebisuya and Briscoe, 2018;Keyte and Smith, 2014).Most of these characterised cases showed that the mechanism responsible for the heterochrony acts downstream of its underlying genetic cause, such as changes in proliferating rates, differences in the initial size of the primordium or distinct rates of protein stability and biochemistry (Gomez et al., 2008;Kicheva et al., 2014;Matsuda et al., 2019;Rayon et al., 2020).Heterochronic shifts can also occur as direct consequence of the causative genetic change, such as those that affect regulatory regions altering the timing of gene expression (Ramaekers et al., 2019).Although differences in gene expression of single transcription factors have the potential to completely modify the subsequent GRN, the relative contribution of such direct heterochrony is in generating morphological diversity remains unknown.Our data indicate that otd is generally expressed more highly during ommatidial maturation in D. mauritiana than D. simulans.
Further analysis shows that otd is actually expressed earlier during ommatidial maturation in D. mauritiana compared to D. simulans.This suggests that cis-regulatory changes in otd lead to ommatidial cells being exposed to Otd for longer in D. mauritiana resulting in larger ommatidia.Together with Ramaerkers and colleagues (Ramaekers et al., 2019), our study shows how morphological diversity in closely related species can be achieved by subtlety altering the temporal expression of a single TF.Importantly, in both cases, these transcription factors, Ey and Otd, act upstream in the GRN controlling the process, thus changes in their expression may promote major differences in downstream effectors.
Further exploration and comparison of the regulatory landscape of otd between D. mauritiana than D. simulans allowed us to identify several candidate cis-regulatory regions that could regulate the eye development of this gene and may contain changes that underlie the differential expression of otd between these two species.These regions may represent novel eye enhancers of otd because they do not overlap with a previously characterised eye enhancer (Hauck et al., 1999;Vandendries et al., 1996), but their activity requires further testing.
We also investigated how these changes in otd expression might alter target gene expression to change ommatidia size.We identified a set of genes that are differentially expressed between these two species when the ommatidia are acquiring their final size that may be acting downstream of Otd, as they have accessible chromatin regions containing putative Otd binding motifs.We compared this set of genes to known and putative targets of Otd which have been characterised later in eye development during pupal stages (Fichelson et al., 2012;Ranade et al., 2008).This comparison showed that a subset of genes for with altered expression in otd mutants are also differentially expressed between D. mauritiana and D. simulans in late L3.In particular, several genes involved in phototransduction (e.g. rh3,slo,Slob,ninaG,inaD,ninaA), genes encoding cytoskeleton and adhesion proteins (Act88F), and other TFs (Dve, vnd, MED10, etc., Suppl.Table 6).
This further suggests that the network downstream of Otd varies between these two species and that ultimately, these changes in the GRN promote differences in ommatidia size between D. mauritiana and D. simulans.Intriguingly, otd also regulates rhodopsin expression in D. melanogaster (Tahayato et al., 2003), and we have shown that several rhodopsins differ in their levels and spatial expression between D. mauritiana and D.
simulans (Hilbrant et al., 2014;Posnien et al., 2012).We showed that Rh3 + ommatidia tend to be larger in diameter than Rh3 -ommatidia.Since rhodopsin 3 expression is higher in D. mauritiana (Hilbrant et al., 2014;Posnien et al., 2012) and Otd directly activates transcription of rhodopsin 3 (Tahayato et al., 2003), a direct link between differences in rhodopsin expression and ommatidia diameter may exist.Therefore, it is possible that differences in otd expression between these two species causes differences in rhodopsin expression as well as ommatidial size, which has important implications for the vision of these flies: increased ommatidium diameter and higher rhodopsin expression could both increase photon capture, potentially resulting in greater contrast sensitivity in eyes with higher otd expression.Conversely, acuity is likely to be reduced due to the inherent tradeoff with sensitivity.This trade-off is heavily influenced by various aspects of visual ecology, such as habitat type, circadian activity patterns, and lifestyle.Thus, substantial functional consequences with strong ecological implications could be linked to changes in the expression of a single gene such as otd.

Conclusions
Our data suggest that changes in the timing of otd expression underlie differences in ommatidia size and thus overall eye size between D. mauritiana and D. simulans.Our work provides new insights into ommatidia size regulation and the evolution of eye size.
Together with evidence from other studies showing that changes in the timing of ey expression contributes to differences in ommatidia number in Drosophila (Ramaekers et al., 2019), we now have a better understand the genetic and developmental mechanisms that underlie the large diversity in Drosophila eye size (Arif et al., 2013;Buchberger et al., 2021;Gaspar et al., 2020;Hilbrant et al., 2014;Keesey et al., 2019;Norry and Gomez, 2017;Posnien et al., 2012;Ramaekers et al., 2019;Reis et al., 2020).Moreover, this evidence suggests that changes in the temporal expression of upstream TFs is a widespread mechanism responsible for morphological evolution.What is also clear is the the potential this system has to build on our existing knowledge of Drosophila eye development (Casares and Almudi, 2016;Casares and McGregor, 2021;Domínguez and Casares, 2005;Gaspar et al., 2019;Kumar, 2018) to ultimately better understand the specification and evolution of organ size more generally.
Female flies of the genotype Otd[YH13], neoFRT19A/FM7c were crossed with males of the genotype RFP, neoFRT19A; ey-Flp.Female F1 progeny were examined for the lack of the Fm7c balancer and these flies were prepared for SEM analysis (see below).

Synchrotron radiation microtomography
Fly heads were removed from the body and placed into fixative (2% PFA, 2.5% GA in 0.1 M sodium cacodylate buffer over night at 4°C.Heads were washed in water, then placed into 1% osmium tetroxide for 48 hours at 4°C, then washed and dehydrated in increasing concentrations of ethanol up to 100%.Heads were then infiltrated with increasing 812 Epon resin concentrations up to 100 % over 5 days and polymerised in embedding moulds for 24 hrs at 70°C. . CC-BY 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is The copyright holder for this preprint this version posted March 17, 2021.; https://doi.org/10.1101/2021.03.17.435774 doi: bioRxiv preprint Heads were scanned at the TOMCAT beamline of the Swiss Light Source (Paul Scherrer Institute, Switzerland; (Stampanoni et al., 2006).Scans were performed using a 16 keV monochromatic beam with a 20 µm LuAG:Ce scintillator.Resin blocks were trimmed and mounted using soft wax and scanned using 20x combined magnification (effective pixel size 325 nm) and a propagation distance of 25 mm.Two thousand projections were taken as the heads rotated through 180°, each with 200 ms exposure.

SEM microscopy
Fly heads were fixed in Bouin's for 2 hours.After 2 hrs, 1/3 of total volume was replaced by 100% ethanol to fully immerse heads in Bouin's and left to fix overnight.Heads were washed and dehydrated 2x 70% EtOH overnight, 2x in 100% ethanol and finally critical point dried and mounted onto sticky carbon tabs on SEM stubs, gold coated and imaged in a Hitachi S-3400N SEM with secondary electrons at 5kV.

Markers and Introgression lines
Males were collected at backcross 7 of three replicate introgression lines (IL1, IL3 and IL4) that were recombinant within the introgressed region (males with phenotypes: yf or vf).
These individuals were genotyped with eleven new additional markers (Suppl.Table 2).
Significant association between each marker and eye size was tested (F-test, type III sum of squares SS) by performing a single-marker ANOVA on the residuals of eye area regressed onto T1 tibia length for each replicate (introgression line (IL 1,3 and 4; n = 20 -60, Suppl.Table 2).Multiple testing was corrected for using Bonferroni correction.All ANOVA models were fitted in the R statistical environment (R Development Core Team 2012) using the CAR package (Fox and Weisberg, 2010).
To narrow down the 2 Mb region the X chromosome region between y and v from D. mauritiana TAM16 into D. simulans YVF was re-introgressed as in (Arif et al., 2013) yf females were backcrossed from multiple replicate lines to yvf males for a further nine generations and the end of the egg-laying cycle of that generation, we collected mothers and genotyped them for molecular markers located in the 2 Mb region (Suppl.Table 2).
Four mothers with breakpoints within this region were identified.Two of them were siblings (IL9.1a and IL9.1b) and they had the same 4 th great-grandmother as IL9.3 and the same .CC-BY 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is The copyright holder for this preprint this version posted March 17, 2021.; https://doi.org/10.1101/2021.03.17.435774 doi: bioRxiv preprint 7 th great-grandmother as IL9.2.Male progeny available for each of these females was collected and genotyped and phenotyped for eye area, ommatidia diameter, ommatidia number and T1 tibia length as described previously in Posnien et al. (Posnien et al., 2012)(Suppl.Table 2).To determine if the D. mauritiana DNA in the 2 Mb region resulted in larger eyes and larger ommatidia, yf males (with some D. mauritiana DNA in the 2Mb interval) were compared to that of their yvf sibling males (i.e., without D. mauritiana DNA) for each introgression line using one-tailed two-sample equal-variance t-tests.

Rh3 related ommatidia measurements
Flies of both species were raised at 25°C with a 12h:12h and dark:light cycle and 40-60% humidity.Density was controlled by transferring 30-40 larvae 20-22 hours after egg laying into fresh food vials.All measurements were performed using female flies 3-5 days after eclosion.
Adult flies were decapitated and heads were cut in half.The halves were fixated in 4% paraformaldehyde (PFA) for 2-3 days.The specimens were washed three times with PBS-T (1xPBS+0.3%Triton X-100) for 10 minutes each.Afterwards they were incubated in 5% H 2 O 2 for 1 day and subsequently in 10% H 2 O 2 for 3-4 days until the heads were depigmented.The depigmented heads were washed three times with PBS (10X: 1.37 M NaCl, 27 mM KCl, 100mM Na 2 HPO 4 , 18mM KH 2 PO 4 ) for 10 minutes each and prepared for immunostaining by incubating them in blocking solution (PBS-T+0.2%goat serum) for 4 hours.The heads were incubated for 48 hrs with the primary antibody (mouse α-Rh3; 1:10 dilution in blocking solution) and then washed five times with PBS for 1 hr each.The secondary antibodies (goat α-mouse Alexa Fluor 555; 1:1000 in blocking solution or goat α-mouse Alexa Flour 647; 1:500 in blocking solution) were applied 48 hours before the heads were washed 3 times with PBS for 10 minutes each.Subsequently, the heads were dehydrated in a graded ethanol series for 30 minutes (30%, 50%, 70%, 90%, 95% and 3 times 100% ethanol in water, respectively) and stored in 100% ethanol.For imaging heads were incubated in methyl-salicylate for 1 hour and then mounted on a cover slip with the eyes facing upwards and imaged with a cLSM (Zeiss LSM 710).The autofluorescence signal was recorded at 488 nm and the staining at 555 nm or 647 nm, respectively.
After acquisition of z-stacks using the cLSM, the data were split into the two channels, converted into tiff-files using Fiji (Image J 1.52n) and analyzed using the 3D image processing program Amira (version 5.4.1).To ensure an accurate threedimensional representation of the data, the voxel size in the z-direction was set to 1.523 to match the refractive index of the cover slips (D 263 M borosilicate glass with refractive .CC-BY 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is The copyright holder for this preprint this version posted March 17, 2021.; https://doi.org/10.1101/2021.03.17.435774 doi: bioRxiv preprint index nD=1.523).Voxel sizes of x and y were entered according to the resolution of acquisition.Volume renderings of the two channels were used to visualize the data and identify ommatidia that were Rh3-positive (Rh3 + ) and those that did not express Rh3 (Rh3 - ) (Figure 4a).Landmarks were used to label the ommatidia and the 3D-measuring tool was used to measure the diameter of the lens of each ommatidium.For landmarks and the measuring tool to work, an isosurface of the autofluorescence channel was created to provide a reference as to where in the 3D space these were placed (Figure 4b).Using this technique, ommatidia diameters of Rh3 + and Rh3 -ommatidia were measured in a central region of the compound eye.Also, in a central region, a quadrant of 10x10 ommatidia was defined and the number Rh3 + and Rh3 -ommatidia was counted (Figure 4b).All raw diameter measurements and ommatidia counts are available in the supplementary material (Suppl.Table 7).

RNA-seq
Flies were raised at 25ºC with a 12h:12h dark:light cycle and their eggs were collected in 2 h time periods.Freshly hatched L1 larvae were transferred into fresh vials in densitycontrolled conditions (30 freshly hatched L1 larvae per vial).Eye-antennal imaginal discs were dissected at three different developmental time points: 72 h after egg laying (AEL), 96 h AEL and 120 h AEL and stored in RNALater (Qiagen, Venlo, Netherlands).Three biological replicates for each sample were generated.Total RNA was isolated using RNeasy Mini Kit (Qiagen).RNA quality was determined using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) microfluidic electrophoresis.
Library preparation for RNA-seq was performed using the TruSeq RNA Sample Preparation Kit (Illumina, catalog ID RS-122-2002) starting from 500 ng of total RNA.Accurate quantitation of cDNA libraries was performed using the QuantiFluor™dsDNA System (Promega, Madison, Wisconsin, USA).The size range of final cDNA libraries was determined by applying the DNA 1000 chip on the Bioanalyzer 2100 from Agilent (280 bp).cDNA libraries were amplified and sequenced using cBot and HiSeq 2000 (Illumina): only 120h eye-antennal imaginal disc samples were sequenced as paired-end (PE) reads (2 x 100 bp), all the rest of samples were sequenced in single-end (SE) reads (1 x 50 bp).
Sequence images were transformed to bcl files using the software BaseCaller (Illumina).
The bcl files were demultiplexed to fastq files with CASAVA (version 1.8.2).The reciprocally re-annotated references described in (Torres-Oliva et al., 2016) were used to map the species-specific reads.Bowtie2 (Langmead and Salzberg, 2012) was used to map the reads to each reference (-very-sensitive-local -N 1) and the idxstats command from SAMtools v0.1.19(Li et al., 2009) was used to summarize the number of mapped reads.HTSFilter (Rau et al., 2013) was used with default parameters to filter out genes with very low expression in all samples.For the remaining genes in each pair-wise comparison, differential expression was calculated using DESeq2 v1.2.7.with default parameters (Love et al., 2014).

ATAC-seq library preparation and sequencing
Samples were obtained following the same procedure as for the RNA-seq experiments: flies were raised at 25 o C with a 12h:12h and dark:light cycle.Freshly hatched L1 larvae were transferred into vials with density-controlled conditions.Eye-antennal imaginal discs were dissected at 96 h AEL and 120 h AEL and maintained in ice cold PBS.Imaginal disc cells were lysed in 50 μl Lysis Buffer (10 mM Tris-HCl, pH = 7.5; 10 mM NaCl; 3 mM MgCl 2 ; 0.1% IGEPAL).Nuclei were collected by centrifugation at 500 g for 5 min.75,000 nuclei were suspended in 50 μl Tagmentation Mix [25 μl Buffer (20 mM Tris-CH 3 COO -, pH = 7.6; 10 mM MgCl 2 ; 20% Dimethylformamide); 2.5 μl Tn5 Transposase;22.5 μl H 2 O] and incubated at 37 C for 30 min.After addition of 3 μl 2 M NaAC, pH = 5.2 DNA was purified using a QIAGEN MinElute Kit.PCR amplification for library preparation was done for 14 cycles with NEBNext High Fidelity Kit; primers were used according to (Buenrostro et al., 2013).Paired end 50 bp sequencing was carried out by the Transcriptome and Genome Analysis Laboratory Goettingen, Germany.

ATACseq peak calling and differential binding site analysis
ATAC-seq raw reads were generated from the following samples (2 replicates each): D.
simulans larvae at 96h AEL and 120h AEL and D. mauritiana larvae at 96h AEL and 120h AEL.These reads were mapped to strain-specific genomes of D. mauritiana and D.
Duplicates were removed using Picard (version 2.20.2) with the parameter REMOVE_DUPLICATE=TRUE.Bam files were then converted to bed files using the Bedtools (version 2.24) bamtobed command.Reads were centred according to (Buenrostro et al., 2013).These reads were then converted to the D. melanogaster coordinate system using liftOver (1.14.0) with custom prepared chain files, one for the conversion of D. mauritiana coordinates to D. melanogaster coordinates and one for the conversion of D. simulans coordinates to D. melanogaster coordinates.Peaks were then called using MACS2 (version 2.1.2,(Zhang et al., 2008)) with the following parameters: -shift -100, extsize 200, -q 0.01.We used the Diffbind package (version 2.12.0, (Ross-Innes et al., 2012)) in R (version 3.6.1.)to search for differentially accessible ATAC-seq regions.A consensus peak set of 19,872 peaks (96h AEL) and 15,868 peaks (120h AEL) was used for all samples and the reads were counted for each identified peak with the dba.count command.For each time point separately we used the dba.analyze command with default parameters to get differentially accessible peaks between the two species.This command uses by default the DESeq2 analysis.All plots were as well generated with the DiffBind package.
To search for TFBM of potential otd regulators, we used the JASPAR core database and its online tool for screening DNA sequences of ATAC-seq peaks with all possible TFB motifs from insects (153 profiles) with a relative profile score threshold of 85%.

Gene regulatory network reconstruction
To search for TFBM of potential otd regulators, we used the JASPAR core database and its online tool for screening DNA sequences of ATAC-seq peaks in the otd locus (between the two flanking genes) with all possible TFB motifs from insects (153 profiles) with a relative profile score threshold of 90%.
To define a list of potential Otd target genes, we used an Otd-motif (Dmelanogaster-FlyFactorSurvey-Oc_Cell_FBgn0004102) from the MotifDB package (version 1.16.1), which provides a collection of available transcription factors in R (version 3.3.3).We searched for Otd binding sites in accesible chromatin regions with the findMotifsGenome.plcommand implemented in the HOMER (version V4.10.4,(Heinz et al., 2010)) in all samples.All peaks with a predicted Otd motif were annotated to an associated gene using the annotatePeaks.plcommand by HOMER and combined all time points and both species into one file.We then looked for the number of genes with an annotated Otd motif and found 1,148 unique genes, which we overlapped with our RNA-seq dataset to find out which of these target genes were differentially expressed between the two species.GO term enrichment analysis of putative Otd target genes was performed using the online tool Metascape (Zhou et al., 2019).We used the online STRING database that integrates all known and predicted associations between proteins based on evidence from a variety of sources (Szklarczyk et al., 2020), to construct networks of DEG encoded proteins.To visualize the network and map genes/prot with Otd motifs we applied the Cytoscape software (Shannon et al., 2003).

In situ hybridisation and immunohistochemistry
In situ hybridizations were carried out using a standard protocol with DIG-labeled antisense RNA probes.Eye-antenna imaginal discs were dissected and fixed at 120 h AEL for 40 min in 4% formaldehyde.To be able to compare the expression patterns avoiding technical differences (i.e.probe affinity and probe concentration), we first aligned the sequences from D. mauritiana and D. simulans and designed RNA probes within fragments with at least 95% of similarity between them (Suppl.Table 8).This design allowed us to perform the in situ hybridization experiments using the same specific probes for each of the candidate genes at the same concentration for both species.The nitro blue tetrazolium/5-bromo-4-chloro-3′-indolyphosphate (NBT-BCIP) reaction was stopped at the same time.Candidate gene sequences were cloned into a TOPO PCR4 (spirit, otd, Ppt1, CG1632, Es2 and CG12112) or pCRII (CG1885, Sptr) vectors (Invitrogen) using specific primer pairs (Suppl.Table 8), respectively, following the manufacturer's protocol.M13 forward and reverse primers were used to linearize the DNA.According to the vector and orientation of the fragments T3, T7 or SP6 RNA polymerase were used to generate the DIG-labeled riboprobes.

Otd positive cells measurements
We used an Analysis of Covariance (ANCOVA) to test for differences between strains for Otd-positive ommatidia while adjusting for differences in development stage by using the number of ommatidial rows as a proxy for the latter.The ANCOVA was performed using base R v4.0.2 (R Core Team, 2020).

Figure 1 .
Figure 1.3D reconstruction and ommatidia size measurements from SRµCT data of female D. simulans (left) and D. mauritiana (right).Facet areas of the ommatidia highlighted in the antero-ventral, central and dorsal-posterior region of the eye are plotted in corresponding colours (far right).Ommatidia number is 996 for the D. simulans and 1018 for the D. mauritiana.Scale bar is 100 µm.

Figure 2 .
Figure 2. Differential and spatial gene expression.(a) Fine-scale mapping of X chromosome QTL.Marker-phenotype association in male recombinant progeny from three replicate introgression lines (IL1, 3, and 4, single-marker ANOVA analysis).Red dashed line indicates the Bonferroni corrected significance threshold of 0.05.Shaded grey area represents a conservative interval of ~ 2 Mb encompassing the X linked QTL.Recombination breakpoints of the new introgression lines (IL9.1-9.3) on the X chromosome (shown for D. simulans Flybase R2.02) define the maximum 662 kb candidate region.White, black, and grey boxes indicate DNA regions from D. simulans YVF, D. mauritiana TAM16 or not determined, respectively (the latter define the maximum candidate region).The table indicates the number of protein coding genes that are present in the candidate region in D. simulans and D. melanogaster.Distribution of eye area (left) and ommatidia diameter (right) measurements by genotype and introgression line.Asterisks indicate significance differences between genotypes where p <.05 (Suppl.Table 2).(b) Differential expression of 49 protein coding genes located in the introgressed region from (a) and expressed at 72 h AEL, 96 h AEL and 120 hAEL.Genes which expression is significantly higher at 120 h AEL in D. simulans are highlighted in blue.Genes significantly upregulated in D. mauritiana at 120 h AEL are shown in red.(c) Expression of differentially expressed genes at 120 h AEL in L3 eye-antenna imaginal discs of D. simulans and D. mauritiana.Open arrowheads indicate the MF, a: antenna, e: eye in (c).

Figure 3 .
Figure 3. otd expression in L3 eye imaginal discs.(a-d) otd mRNA at 110 h AEL and 120 h AEL in D. mauritiana (a-b) and D. simulans (c-d).Arrowhead indicates the MF.Asterisks indicate expression in the ocellar region.D. mauritiana already exhibits otd mRNA at 110h (red arrowhead).(e-f')Immunostaining showing Otd protein (magenta, e') in mature ommatidia (marked in green by Elav) and the ocellar region (marked by an asterisk) in D. mauritiana (e-e') and D. simulans (f-f').Staining against actin (in blue) was used to mark the MF.(g) Plot showing the number of Otdpositive ommatidia rows (X axis) at different developmental time points (y axis, developmental points infer by number of ommatidia rows).Asterisks represent statistical significance p< 0.001.Scale bars: 50 µm.
in the otd locus during eye development between D. simulans and D. mauritiana Our mapping and expression analyses indicate that the differences in otd expression likely contribute to differences in ommatidia size between D. simulans and D. mauritiana.Given that there is only 1 amino acid difference in the Otd sequence between our focal strains of D. simulans and D. mauritiana, our data suggest that the causative changes are located in otd regulatory regions.Due to the microsyntenic conservation between D. melanogaster, D. simulans and D. mauritiana, we considered the regulatory landscape of otd as the region between its two flanking genes, Caf1-180 and CG12772, revealed by the presence of a Topological Associated Domain (TAD) in the corresponding D. melanogaster region (a region of 69 kb in D. mauritiana and 70 kb D. simulans, Fig. 4a, Suppl.Fig. 4; http://chorogenome.ie-freiburg.mpg.de/).

Figure 4 .
Figure 4. Chromatin accessibility at the otd locus.(a) Open chromatin peaks at the otd locus in 96 h AEL and 120 h AEL D. mauritiana and D. simulans eye-antenna imaginal discs.(b) Detail of differential peaks: 6, 11, 16, 17 and 19.(c) Alignment of the sequence of the first D. mauritiana specific region in peak 16 with D. simulans.(d) Alignment of the sequence of the second D. mauritiana specific region in peak 16 with D. simulans.
antennal disc development.Peak 11 encompasses 412 putative TFBMs in D. mauritiana and 411 putative sites in D. simulans, including Mad, So, Optix, Exd motifs exclusively in D. mauritiana, and Trl and Ttk predicted specifically in D. simulans.For the D. mauritianaspecific peak 16, 37 TFBM were predicted in the first D. mauritiana-specific peak (peak 16a, threshold 85%) for TFs such as Brk, Caup and Mirror, Clamp and Trl but only 8 TFBMs were predicted for the second D. mauritiana-specific peak (peak 16b, threshold 80%) including sites for CTCF, Mad, Brk and Trl.Peak 17 contained 269 predicted TFBMs in D. mauritiana and 275 predicted TFBMs in D. simulans, including Dichaete (D) and CTCF motifs, only present in the later species.Regarding peak 19, which covers the longest sequence, 790 and 754 TFBM were predicted for D. mauritiana and D. simulans, respectively.From those, 57 were D. mauritiana-specific (Eip74EF, Inv, Dll, Mirr, Al, Lim1, Vvl, Suppl.Table

Figure 5 .
Figure 5. Otd downstream targets.(a) 161 genes upregulated in D. mauritiana have an associated peak that contains at least one Otd motif.(b) GO enrichment for those genes with Otd motifs in D. mauritiana (c) 111 genes upregulated in D. simulans have an associated peak that contains Otd motif.(d) GO enrichment for those genes with Otd motifs in D. simulans.

Figure 6 .
Figure 6.Differences in size of Rh3 positive ommatidia.(a) Three-dimensional reconstruction (volume rendering) of a representative D. simulans compound eye after clearing and cLSM imaging showing Rh3 expressing ommatidia.a -anterior, p -posterior, d -dorsal, v -ventral.(b) Isosurface reconstruction based on the autofluorescence (volume rendering).The 10x10 central ommatidia that were counted are labelled with yellow dots.(c) Boxplot showing the diameter of Rh3 + and Rh3 -ommatidia in both species.Differences in diameter were statistically tested applying one-way ANOVA (F 3,362 = 8.571, P = 0.0000164) followed by Tukey multiple comparisons (Pvalues).(d) Boxplot showing the percentage of Rh3 + ommatidia based on the 10x10 central ommatidia counted.The differences are not significantly different (χ 2 = 3.0523, df = 1, p-value = 0.08062).The raw data for c and d are available in Suppl.Table7.
was used to drive expression of the transgenes.To generate mitotic clones of mutant otd in developing eyes we used the stocks Otd[YH13], neoFRT19A/FM7c and RFP, neoFRT19A; ey-Flp which were obtained from Bloomington Stock Centre (Stock nos.#8675 and #67173 respectively). .
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is . CC-BY 4.0 International license made available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It isThe copyright holder for this preprint this version posted March 17, 2021.; https://doi.org/10.1101/2021.03.17.435774 doi: bioRxiv preprint FastQC software (version 0.10.1,Babraham Bioinformatics) was performed.All RNAseq reads are accessible in the Short Read (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is Before the mapping step, PE 100 bp reads were converted into SE 50 bp by splitting the reads in half and merging right and left reads into a single file.
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is