Gene protein sequence evolution can predict the rapid divergence of ovariole numbers in Drosophila

Ovaries play key roles in fitness and evolution: they are essential female reproductive structures that develop and house the eggs in sexually reproducing animals. In Drosophila, the mature ovary contains multiple tubular egg-producing structures known as ovarioles. Ovarioles arise from somatic cellular structures in the larval ovary called terminal filaments, formed by terminal filament cells and subsequently enclosed by sheath cells. As in many other insects, ovariole number per female varies extensively in Drosophila. At present however, there is a striking gap of information on genetic mechanisms and evolutionary forces that shape the well-documented rapid interspecies divergence of ovariole numbers. To address this gap, here we studied genes associated with D. melanogaster ovariole number or functions based on recent experimental and transcriptional datasets from larval ovaries, including terminal filaments and sheath cells, and assessed their rates and patterns of molecular evolution in five closely related species of the melanogaster subgroup that exhibit species-specific differences in ovariole numbers. From comprehensive analyses of protein sequence evolution (dN/dS), branch-site positive selection, expression specificity (tau) and phylogenetic regressions (PGLS), we report evidence of 42 genes that showed signs of playing roles in the genetic basis of interspecies evolutionary change of Drosophila ovariole number. These included the signalling genes upd2 and Ilp5 and extracellular matrix genes vkg and Col4a1, whose dN/dS predicted ovariole numbers among species. Together, we propose a model whereby a set of ovariole-involved gene proteins have an enhanced evolvability, including adaptive evolution, facilitating rapid shifts in ovariole number among Drosophila species. Significance Statement Ovaries in Drosophila, like in other insects, contain egg producing structures, known as ovarioles. The number of ovarioles per female varies among Drosophila species, but little is known about the genes and evolutionary dynamics that may shape interspecies changes in ovariole numbers. Here, used a priori experimental and transcriptome data from D. melanogaster to identify genes involved in ovariole formation and functions, and studied their molecular evolution among its closely related species within the melanogaster subgroup. Using a multi-layered analysis consisting of protein sequence divergence (dN/dS), adaptive evolution, expression breadth, and phylogenetic regressions, we identified 42 genes whose molecular evolution patterns were well linked to ovariole numbers divergence. Further, gene protein sequence divergence was often predictive of species ovariole numbers.


22
Ovarian development is a process that is poised to play key roles in organismal evolutionary 23 biology, as the female gonads form and house the oocytes and/or eggs that are central to fertility and 24 reproductive success of a species, and thus affect their fitness (Miller et al. 2014;Macagno et al. 2015). 25 In insects, the most well-studied model with respect to ovarian development and genetics is the fruit fly 26 Drosophila melanogaster (Dansereau and Lasko 2008;Eliazer and Buszczak 2011;Li et al. 2014;Slaidina 27 et al. 2020;Lebo and McCall 2021). The mature ovary in D. melanogaster, as in other species of insects, 28 is comprised of tubular egg-producing structures known as ovarioles (King et al. 1968;Dansereau and 29 Lasko 2008; Lebo and McCall 2021), which are a central factor shaping organismal reproductive output 30 (Montague et al. 1981;Starmer et al. 2003;Church et al. 2021). The number of ovarioles contained in the 31 ovaries is highly variable within the genus Drosophila (Kambysellis and Heed 1971;Hodin and Riddiford 32 2000; Starmer et al. 2003;Markow et al. 2009;Sarikaya et al. 2019;Church et al. 2021). As an example, 33 within the melanogaster subgroup, D. melanogaster has typically about 19 ovarioles per ovary, while its 34 closely related sister species D. sechellia has only about 8 to 9 ovarioles per ovary (Hodin and Riddiford 35 2000). A broad range of ovariole numbers has been observed across the family Drosophilidae, from one 36 to more than 50 per ovary across the genus Drosophila (Sarikaya et al. 2019;Church et al. 2021). At 37 present, however, we know little about the genetic basis of the evolution of ovariole number within insects 38 (Hodin and Riddiford 2000;Markow et al. 2009;Sarikaya et al. 2019). 39 A central factor that may underlie the rapid interspecies transitions in ovariole numbers in 40 Drosophila is the evolvability of ovariole-related protein-coding genes, that is, the propensity of the 41 proteins encoded by these genes to diverge and/or undergo adaptive sequence changes (Wagner and Zhang 42 2011; Cutter and Bundus 2020). Functional amino acid changes in protein-coding DNA and associated 43 selection pressures (measured as nonsynonymous to synonymous changes, or dN/dS (Yang 1997;44 to posterior direction between the TFs, depositing basement membrane that partitions the remaining cells 91 of the ovary (germ cells and posterior somatic cells) into the developing ovarioles (King et al. 1968;King 92 1970; Slaidina et al. 2020). Additional somatic cells in the LL3 ovary include intermingled cells, which 93 are interspersed between the germ cells and are involved in their proliferation (Gilboa and Lehmann 2006), 94 cap cells, which form the adult germ line stem cell niche (Song et al. 2002), follicle stem cell precursors, 95 which give rise to adult follicle stem cells (Slaidina et al. 2020;Slaidina et al. 2021), and swarm cells, 96 whose precise functions largely remain to be ascertained (Slaidina et al. 2020) (fig. 1A). In this regard, 97 understanding the interspecies evolution of ovariole number in Drosophila requires consideration of the 98 genes and proteins regulating cell behaviour in the larval ovary, and particularly the behaviours of the TF 99 and SH cells, which are instrumental to determining ovariole numbers in D. melanogaster. 100 Until recently, research on the relationships between divergence in gene sequences and ovariole 101 numbers in Drosophila was challenged by the lack of data on the identity of protein-coding genes 102 expressed in somatic cells of the larval ovary that regulate ovariole number (Sarikaya et al. 2012;Sarikaya 103 and Extavour 2015). Recently available large-scale functional genetic and cell type-specific expression 104 data from D. melanogaster, however, now provide a means to systematically identify genes linked to 105 ovariole numbers, and in turn, assess their molecular evolution across species. A large-scale RNAi screen 106 of 463 signalling genes from 14 conserved animal signalling pathways revealed that TF-mediated ovariole 107 number determination is regulated by all conserved animal signalling pathways (Kumar et al. 2020). 108 Another study using bulk-RNA seq expression data from FACS-separated germ cells and somatic cells 109 revealed additional genes differentially expressed throughout TF formation, suggesting their potential 110 involvement in ovariole number regulation (Tarikere et al. 2022). In addition to those studies, a recent sc-111 RNA seq study yielded unique transcriptional profiles for all of the known cell types in the D. 112 melanogaster LL3 ovaries ( fig. 1), providing a novel resource to identify and study the evolution of genes 113 transcribed in TF and SH cells, the two crucial cell types in determining ovariole number (Slaidina et al. 114 2020). 115 Collectively these datasets provide valuable empirical data from which to a priori identify sets of 116 genes that regulate ovariole numbers or functions in Drosophila, and in turn, to evaluate which of these 117 genes exhibit elevated or otherwise unusual rates of interspecies protein sequence evolution, including 118 adaptive evolution, suggesting them as candidates for driving interspecies divergence of ovariole numbers 119 in Drosophila. For example, by assessing dN/dS, we may ask whether ovariole-related gene protein 120 sequences typically have been under strict purifying selection, which could mean that phenotypes 121 regulated by these genes are likely to show high pleiotropy and low evolvability, and to have minimal 122 potential to diverge neutrally or adaptively (Fisher 1930;Otto 2004;Wagner and Zhang 2011;Cutter and 123 Bundus 2020; Munds et al. 2021). If, in contrast, some ovariole-related genes have been subjected to 124 relaxed selection and/or have commonly experienced adaptive changes, we might expect high phenotypic 125 evolvability and adaptability (Otto 2004;Larracuente et al. 2008;Clark et al. 2009;Mank and Ellegren 126 2009; Montgomery et al. 2011;Luke et al. 2014;Corso et al. 2016;Chebbo et al. 2021). In this regard, 127 the study of the evolution of protein-coding genes (from dN/dS) that are pre-screened for likely roles in 128 ovariole numbers and/or functions by studies like the ones described above (Kumar et al. 2020;Slaidina 129 et al. 2020;Tarikere et al. 2022) provides a novel pathway to advance our understanding of the genetic 130 factors and evolutionary forces that shape rapid interspecies divergence in ovariole numbers. 131 In the present study, we rigorously assess the molecular evolutionary patterns of genes that regulate 132 ovariole numbers and/or functions, identified based on one or both of functional genetic evidence (Kumar 133 et al. 2020) or transcriptional activity (Slaidina et al. 2020;Tarikere et al. 2022). We focus on the 134 molecular evolution of ovariole-related genes within five species of the melanogaster subgroup of 135 Drosophila, a closely related species clade that includes D. melanogaster, diverged from a common 136 ancestor about 12.6 mya (Tamura et al. 2004), and exhibits substantial interspecies variation in ovariole 137 numbers (Hodin and Riddiford 2000;Starmer et al. 2003;Markow et al. 2009). We identify 42 genes that 138 are high confidence candidates for contributing to the genetic basis of interspecies divergence in ovariole 139 numbers. We hypothesize that evolved changes in these genes are apt to underlie ovariole number 140 divergence among taxa given that they exhibit an ovariole-related function (Kumar et al. 2020;Slaidina 141 et al. 2020;Tarikere et al. 2022), have a propensity to diverge in protein sequence, or high evolvability, 142 show a high frequency of adaptive sequence evolution events in branches of the phylogeny, and are 143 associated with low pleiotropy (Yanai et al. 2005). Further, phylogenetic regressions show gene dN/dS 144 has predictive associations to ovariole numbers. Collectively, our findings provide a genetic framework 145 to explain the rapid interspecies divergence of ovariole numbers in Drosophila, which we propose is 146 largely mediated by selection pressures shaping the evolution of functional protein sequences, and thus 147 ovariole numbers. group was used as an outgroup for phylogeny construction, see "Drosophila Phylogeny" section; the 155 abbreviated names were used in tables and figures). The melanogaster subgroup had the following 156 advantages for our study: (1) each species has a well-annotated whole genome sequence available for D. simulans (33.9), D. yakuba (25.8) and D. erecta (27.0) ( fig. 2; see values and mild variability (Hodin 160 and Riddiford 2000; Starmer et al. 2003;Markow et al. 2009); (3) the close relatedness of the five species 161 (Tamura et al. 2004;Cutter 2008) minimizes biological differences other than ovariole numbers among 162 taxa, a feature that facilitates detection of cause-effect relationships (here, dN/dS and ovariole number 163 (Felsenstein 1985;Bromham et al. 1996;Whittle and Johnston 2003;Thomas et al. 2010); (4) the taxa 164 have known high confidence orthologs (Waterhouse et al. 2011;Gramates et al. 2022) and previously 165 calculated M0 dN/dS values (M0 is a single value per phylogeny) per gene from PAML (Yang 1997, 166 2007; Stanley and Kulathinal 2016), as well as five-species codon alignments per gene available for 167 customized dN/dS analysis (Stanley and Kulathinal 2016); (5) the dN and dS values among the species in 168 this subgroup have substantially diverged, yet are also unsaturated in the frequency of substitutions, and 169 thus are within the ideal range for dN/dS analysis (Castillo-Davis et al. 2004;Larracuente et al. 2008;170 Treangen and Rocha 2011) (for example, we found that the 95 th percentile for M0 dN=0.235 and M0 171 dS=0.791 for the 9,232 genes that had orthologs in all five species and M0 values) and; (6) all species in 172 the clade are very closely related to D. melanogaster, the species for which experimental and 173 transcriptome data on genes associated with ovariole roles are available (Kumar et al. 2020;Slaidina et 174 al. 2020;Tarikere et al. 2022). In sum, this taxonomic group is especially well suited to the study of the 175 evolution of ovariole-related genes.

176
The dN/dS value reflects the rate of protein divergence and the potential types of selective 177 pressures that may have affected a gene (Yang 1997(Yang , 2007. Values of dN/dS <1 suggest a history of 178 purifying selection, =1 infer neutral evolution, and >1 suggest a history of positive selection (Yang 1997, 179 2007). However, even when dN/dS <1 across an entire gene (which is a conservative measure of dN/dS 180 (Yang 2007)), elevated dN/dS values in one gene relative to another suggest an enhanced degree of 181 positive selection and/or neutral evolution (Yang 1998(Yang , 2007Buschiazzo et al. 2012;Ho and Smith 2016;182 Mitterboeck et al. 2017;Whittle et al. 2021). The gene-wide dN/dS may be determined as a single M0 183 dN/dS value per clade (Yang 2007) and/or as free-ratios dN/dS with separate values for each branch in a 184 clade (Yang 2007). Further, as gene-wide dN/dS analyses may be innately conservative, a fine-scale 185 branch-site analysis may also be used to test for positive selection at gene codon sites within specific 186 species branches of interest within a tree (Zhang et al. 2005;Yang 2007 ovariole or egg laying phenotypes (and passed screening of z>1; note that genes may belong to more than 199 one category) (Kumar et al. 2020). For these four genes sets, we identified any genes with M0 dN/dS ≥1.5 200 higher than the genome-wide median. The cutoff was marginally lower than the other two datasets 201 described below because of the innate conserved nature of these signalling pathway genes, which are 202 largely at least as old as animal divergence, in excess of 600 million years (Srivastava et al. 2010;Kumar 203 et al. 2020) (see additional details in Supplementary Text File 1).

204
The second is the BULKSG dataset, based on bulk-RNA seq data obtained from pooled larval 205 ovarian somatic cells or germ cells from the early (72 hours after egg laying = 72h AEL), mid (96h AEL) 206 and late (120h AEL) TF developmental stages (Tarikere et al. 2022); P-values were from DeSeq2 (Love 207 et al. 2014). For this gene set, we screened for any differentially expressed genes that had M0 dN/dS≥0.20 208 in the melanogaster subgroup for further study. This represents a value ≥2.2 higher than the genome-wide 209 median. For this dataset we chose a higher cutoff than for the SIGNALC dataset, since BULKSG were 210 considering genome-wide transcript data and were not limited to the members of highly conserved 211 signalling pathways. some genes were upregulated in more than one cell type using the criteria therein (Slaidina et al. 2020).

216
For genes with differential expression in one cell type relative to the others (P<0.05), we identified those   For the ovariole-related genes identified as rapidly evolving (with elevated M0 dN/dS) from the 228 SIGNALC, BULKSG, and SINGLEC datasets, and for all genes in the genome, we used the available 229 five-species alignments that had an ortholog in all species in the melanogaster subgroup (Stanley and 230 Kulathinal 2016) and determined the M1 free ratios dN/dS per species terminal branch using codeml in 231 PAML (Yang 2007), which allows a separate dN/dS value for each branch. Codeml is based on maximum 232 likelihood in deriving estimates of dN/dS values, and default parameters were used in the assessments 233 (Yang 2007). The M1 model (Yang 2007) has been commonly and effectively used to study branch dN/dS 234 ( fig. 2) (e.g., (Dorus et al. 2004;Nadeau et al. 2007;Clark et al. 2009;Wlasiuk and Nachman 2010;235 Mensch et al. 2013;Borges et al. 2019;Kong et al. 2019;LaBella et al. 2021)). Using the dN/dS for each 236 of the five terminal species branches, we assessed associations with respect to species transitions in 237 ovariole numbers (terminal species branch analysis), an approach that has proven effective for determining 238 relationships between dN/dS values and phenotypes of interest (Dorus et al. 2004;Nadeau et al. 2007; 239 Wlasiuk and Nachman 2010). respectively had values below this threshold, and we found even higher percentages (up to 100%) for the 245 four other species. Only gene branches that had dN or dS >0.001 were included for further assessment to 246 ensure sufficient divergence for study (Cusack and Wolfe 2007;Whittle et al. 2021). The minority of cases 247 of a branch where dN was >0.001 and dS was at or near zero were denoted simply as "dN/dS>1" (e.g., selection at one or more codon sites in that species branch. We studied the presence or absence of branch-267 site positive selection in each gene, suggested by Zhang et al. (2005), without including the post-hoc 268 option for BEB probability analysis per codon site that has low power (Zhang et al. 2005), which was 269 most informative for our objective of identifying genes that have sites that experienced positive selection. here (Yanai et al. 2005). For this, we accessed expression data from 59 tissue types and developmental 287 stages from D. melanogaster (30 developmental stages and 29 tissues, table S1). The data include gene 288 expression levels (RPKM) across development for embryos (12 stages), larvae (6 stages), pupae (6 stages) 289 and adults (3 stages of males/females), and for major tissue types of the adult males and females (including where n=number of tissues/stages studied, i= tissue/stage, = expression level of gene in tissue/stage i, 299 and max (xi)= the expression level in the tissue/stage type with maximum expression (Yanai et al. 2005).

322
We focused on dN/dS and branch-site selection tests to assess all five species of the melanogaster 323 subgroup for positive selection (Yang 2007). In addition, as a supplemental analysis, we conducted

357
We first report the results for the ovariole-related genes from the SIGNALC dataset (Kumar et al. 358 2020). We found that the ovariole-related SIGNALC genes exhibited very low M0 dN/dS for signalling 359 genes that affected ovariole number and/or egg laying (MWU-tests had P<0.05 versus the genome-wide 360 values; fig. 3A). This suggests a history of strong purifying selection on these highly conserved signalling 361 genes, which may be partly due to their high pleiotropy, given that all of these signalling pathways play 362 multiple roles in development and homeostasis (Kumar et al. 2020 (table 1), which is consistent with potential high adaptability of these genes. Of particular note is the D. 389 sechellia branch, as this species evolved a very low ovariole number (17 ovarioles per female, fig. 2), only 390 half that of its most closely related sister species D. simulans (33.9 ovarioles per female, fig. 2), since 391 diverging from their last common ancestor only about 3 mya (Cutter 2008   The patterns in table 1 support the hypothesis that protein sequence changes, including adaptive 416 changes, in these ovariole-related genes may underlie the genetic basis for the marked divergence in

427
We identified genes whose high differential expression in the D. melanogaster larval ovary 428 suggested a role in ovariole number regulation using the BULKSG RNA-seq datasets using pooled larval 429 ovarian somatic versus pooled germ cells from different stages of TF formation (Tarikere et al. 2022

467
The SINGLEC dataset was based on sc-RNA seq data generated from the late third instar D.  Col4a1 play roles in basement membrane formation (Yasothornsrikul et al. 1997;Kiss et al. 2019), and 491 that SHm cells lay the membrane that separates the TFs for ovariole development (King 1970;Slaidina et 492 al. 2020). Given the crucial roles of these cell types in determining ovariole number (King, Aggarwal et 493 al. 1968, Sarikaya andExtavour 2015), the rapid evolution of these five genes may partially underlie 494 ovariole number divergence between species (King et al. 1968;Sarikaya and Extavour 2015) in the 495 melanogaster subgroup (table 4).

496
In terms of molecular evolution per terminal species branch, the five genes in  The genes identified above as highly expressed in TF and SH cells, could also be highly expressed 527 in additional cell types (Slaidina et al. 2020). Indeed,on average we found that differentially expressed 528 genes were upregulated in 1.9±0.02 cell types. Thus, for additional stringency we isolated the subset of

Functional predictions for upregulated TF and SH genes 548
The studied molecular evolutionary parameters for all genes studied in fig. 4  identify genes putatively involved in the evolution of phenotypes, as has been also observed for other 575 diverse traits across multiple taxa (Dorus et al. 2004;Nadeau et al. 2007;Ramm et al. 2008;Wlasiuk and 576 Nachman 2010; Luke et al. 2014;Corso et al. 2016;Chebbo et al. 2021).  (table 1, table 3, table 4, table S6, fig. 4, fig. S4). This 594 hypothesis is further supported by the fact that all of the ovariole-related genes revealed herein (tables 1-595 5) have been explicitly demonstrated to regulate ovariole number (Kumar, Blondel et al. 2020), and/or are 596 highly and/or exclusively expressed in somatic ovarian cells whose behaviour determines ovariole number 597 (table 1, table 3, table 4, table S6, fig. 4) (King et al. 1968;King 1970;Sarikaya et al. 2012;Sarikaya et al. 598 2019; Slaidina et al. 2020;Tarikere et al. 2022).  , table 1, table 3, table 4). Thus, low pleiotropy may have partly 613 contributed to high evolvability, and enhanced adaptive potential. These events of positive selection in the 614 ovariole-related genes (table 1, table 3, table 4, fig. 4), may have arisen by natural selection for adaption 615 to changes in environment or oviposition substrates (Jagadeeshan and Singh 2007), and/or may have often 616 been driven by the widely-reported and dynamic sexual behaviors of Drosophila, as described below.

619
Sexual selection may contribute to the adaptive evolution of reproductive characteristics and genes 620 in animals (Swanson and Vacquier 2002;Clark et al. 2009), including in Drosophila (Civetta and Singh 621 1998; Swanson et al. 2004;Proschel et al. 2006). Thus, one possibility is that this phenomenon may shape 622 the evolution of ovariole-related genes observed herein (table 1, table 3, table 4 (across species table 1, table 3, table 4), nor to give rise to the observed 658 predictive relationships between dN/dS and ovariole numbers using PGLS (table 5) table 1, table 3, table 4, table S8, fig. 4, fig. S4), is unlikely to be explained by neutral evolution alone.

662
Thus, the present data suggest that neutral evolution has not been the only or main driving factor shaping 663 amino acid changes in ovariole-related genes in the melanogaster group, which we propose instead are 664 best explained by a history of adaptive evolution.

665
Another factor in addition to narrow expression breadth (a factor that affects individual genes) that 666 could in theory lead to relaxed purifying selection on nonsynonymous mutations in ovariole genes is small 667 population size, which may affect entire genomes (Kimura 1962;Strasburg et al. 2011;Gossmann et al.  (table 1, table 2, table 4,table S6 proliferation at a specific rate and to a specific degree during larval stages, morphogenetic movements 694 including intercalation and migration to establish terminal filaments, and extracellular matrix deposition 695 to separate ovarioles from each other within the gonad (King 1970). Any of these developmental processes 696 could in principle be the target of evolutionary change in interspecies ovariole number divergence. Indeed, 697 we previously showed that evolution of different developmental mechanisms underlies convergent 698 evolution of similar ovariole numbers between or within species . Accordingly, 699 we would expect that the genes underlying these evolutionary changes might play roles in multiple 700 different developmental processes, and this prediction is supported by our findings herein. The genes that 701 we have identified here as not only rapidly evolving in this subgroup (table1 , table 3, table 4), but also 702 with molecular evolutionary rates that are highly predictive of lineage-specific ovariole numbers ( subgroup members and to its hypothesized last common ancestor  Our results demonstrate the utility of comprehensive dN/dS and positive selection analyses in 718 identification of rapidly evolving genes that may shape the evolution of a core reproductive phenotype.

719
This method provides valuable opportunities for the discovery of genes and evolutionary processes 720 involved in interspecies phenotype divergence (Dorus et al. 2004;Nadeau et al. 2007;Ramm et al. 2008; in each species for the rapidly evolving ovariole-related genes identified here (table 1, table 3, table 4,   732   table S6)  Hawaiian Drosophila, an expansive taxonomic group that exhibits marked phenotypic diversity in sexual 751 characteristics ranging from behaviours to ovariole numbers (Carson 1997;Singh and Singh 2014), 752 (Sarikaya et al. 2019). Studies on the relationships between protein sequence changes and ovariole 753 numbers in Hawaiian Drosophila will be facilitated by increased collection of genomic data from a range 754 of species in this taxon group, and transcriptomic data for the larval ovaries, including TFs and SH cells.

755
Such research will help further decipher the genetic factors shaping the rapid evolution of ovariole 756 numbers in the Drosophila genus, and thus in insects more broadly.    Kumar et al. (2020) and their expression status in the soma and germ cells in the larval ovary (each pooled across stages), and among somatic cells at the early, mid and late stages (DeSeq2 P<0.01 (Tarikere et al. 2022)). Note that if a gene is designated as upregulated in the germ cells, this automatically indicates it is downregulated in soma. A total of 25 of 27 genes showed differential expression using at least one of these comparisons.

Somatic cellsstage FBgn0011274
Dif Soma Early FBgn0014020 Rho1 Pten Tefu Germ -Notes: the gene Paris (FBgn0031610), that was rapidly evolving in Dmel-Dsim but lacked high confidence orthologs in all five species was reported as upregulated in the germ cells relative to the soma cells in (Tarikere et al. 2022).   (M0 dN/dS>0.20) and that were highly upregulated at one stage of the larval ovary somatic cells (versus the others; three stages early, mid, late, among the top 30 most upregulated genes, table S5) in Dmel using BULKSG data (Tarikere et al. 2022) and that also exhibited upregulation in at least one cell type (versus all others) using SINGLEC data among the nine studied LL3 ovary cell types (Slaidina et al. 2020). Shown are the dN/dS per species branch, the presence of branch-site positive selection (P<0.05), the tau values and an example of key functionality as described in DAVID (Huang da et al. 2009). "Stage up" indicates the larval ovary stage where the gene was upregulated (P<0.05). SingleC up indicates the cell type(s) with upregulation.  Table 5. PGLS analysis of the relationship between ovariole number and dN/dS for genes putatively involved in ovariole number evolution from tables 1, 3 and 4 (42 genes total). The 17 genes that showed a relationship using PGLS are shown (P<0.05), and includes the intercept, the slope, and the predicted ovariole numbers using the model. In addition, the dataset that each gene was identified from and the  fig. 2 and branch lengths used for PGLS analyses are in Materials and Methods.   (Tamura et al. 2021) and DNA sequence data from DrosoPhyla (Finet et al. 2021). The five species of the melanogaster subgroup are shown. The relatively distantly related D. ananassae (Dana) was used as an outgroup for tree construction. Ovariole numbers (ON) are shown and are for two ovaries per female and are from the following sources: D. melanogaster (Dmel), D. sechellia (Dsec), and D. yakuba (Dyak) (Hodin and Riddiford 2000), D. simulans (Dsim) (averaged, (Hodin and Riddiford 2000;Starmer et al. 2003) and D. erecta (Dere) (Markow et al. 2009) (see respective articles for variation). All nodes had 100/100 bootstrap support. Figure. 3. Box plots of A) M0 dN/dS of genes with five-species orthologs in the melanogaster subgroup for each of four groups of signalling/connector genes that affected ovariole/egg numbers using RNAi in D. melanogaster (Kumar et al. 2020) and for the genome-wide values; and B) tau for all genes in each of the four groups of ovariole number/egg laying affecting genes and the genome-wide values. Different letters (a.b) below bars indicate a statistically significant difference (MWU-tests P<0.05) between the genome-wide values and each group of genes. The median and 25 th percentiles are shown for dN/dS and tau as reference points for the genome-wide values (that is, across all 9,232 genes with known dN/dS and five-species orthologs). . Swarm cells (SW) cells were excluded as too few genes were rapidly evolving for study (SW: 4). Note that a gene could be upregulated in more than one cell type. The genome-wide values are for all genes with five-species orthologs in the melanogaster subgroup.