Gene expression differentiation in the reproductive tissues of Drosophila willistoni subspecies and their hybrids

Early lineage diversification is central to understand what mutational events drive species divergence. Particularly, gene misregulation in interspecific hybrids can inform about what genes and pathways underlie hybrid dysfunction. In Drosophila hybrids, how regulatory evolution impacts different reproductive tissues remains understudied. Here, we generate a new genome assembly and annotation in Drosophila willistoni and analyse the patterns of transcriptome divergence between two allopatrically evolved D. willistoni subspecies, their male sterile and female fertile hybrid progeny across testis, male accessory gland, and ovary. Patterns of transcriptome divergence and modes of regulatory evolution were tissue‐specific. Despite no indication for cell‐type differences in hybrid testis, this tissue exhibited the largest magnitude of expression differentiation between subspecies and between parentals and hybrids. No evidence for anomalous dosage compensation in hybrid male tissues was detected nor was a differential role for the neo‐ and the ancestral arms of the D. willistoni X chromosome. Compared to the autosomes, the X chromosome appeared enriched for transgressively expressed genes in testis despite being the least differentiated in expression between subspecies. Evidence for fine genome clustering of transgressively expressed genes suggests a role of chromatin structure on hybrid gene misregulation. Lastly, transgressively expressed genes in the testis of the sterile male progeny were enriched for GO terms not typically associated with sperm function, instead hinting at anomalous development of the reproductive tissue. Our thorough tissue‐level portrait of transcriptome differentiation between recently diverged D. willistoni subspecies and their hybrids provides a more nuanced view of early regulatory changes during speciation.


| INTRODUC TI ON
A central goal in evolutionary biology is finding the genetic alterations underlying ecological, physiological, and morphological change (Lewontin, 1974). Mutations affecting regulation of gene expression can impact how the species blueprint becomes implemented molecularly during life cycle, across cell types, or under particular external cues, mediating phenotypic change during adaptation and species divergence (Glaser-Schmitt & Parsch, 2018;Hagen et al., 2019;Kratochwil et al., 2018). Early stages during lineage differentiation merit special attention as it is when the differential presence of genetic regulatory variants between two decreasingly connected gene pools can influence the consolidation of phenotypic differentiation and speciation. This differential presence of regulatory variants is ultimately reflected in diverging expression, resulting in some cases in disrupted regulatory networks, hybrid misexpression and dysfunction (Civetta, 2016;Cutter & Bundus, 2020;Mack & Nachman, 2017;Signor & Nuzhdin, 2018).
Different methodologies have been used to uncover the regulatory changes underlying changes in mRNA levels (Gilad et al., 2008;Wittkopp et al., 2004). Particularly, allele-specific expression (ASE) studies have revealed pervasive trans-regulatory variation within species, a more preeminent role of cis-regulatory changes in the interspecific evolution of mRNA levels, and normally larger size effects on expression differences in the case of cis-as opposed to trans-acting changes (Coolon et al., 2014;Emerson et al., 2010;Kopania et al., 2022;Zhang et al., 2011). Notably, the interpretation of results from ASE studies involving different species can be distorted by developmental defects that result in cellular tissue composition differences (Brawand et al., 2011;Good et al., 2010;Hunnicutt et al., 2022). Therefore, observations drawn from ASE studies using early divergent lineages in which cellular tissue composition differences are minimal are particularly informative (Mugal et al., 2020).
A key aspect of expression divergence is how lineage differentiation is distributed across the genome, in particular on the sexual chromosomes versus the autosomes. Different studies point towards a larger interspecific differentiation on the X (Brawand et al., 2011;Dean et al., 2015;Llopart, 2012;Meisel et al., 2012b), a pattern known as faster-X evolution, which can result from both neutral and adaptive evolution (Vicoso & Charlesworth, 2009). In Drosophila, adaptive expression evolution underlying the faster-X effect has been reported in embryos (Kayserili et al., 2012), for male-biased genes in adults (Llopart, 2012), and for tissue-restricted expressed genes also in adults (Meisel et al., 2012b). Importantly, the X chromosome contributes disproportionately to postzygotic reproductive isolation between diverging lineages (Coyne, 1992;Haldane, 1922). This prominent role in hybrid dysfunction in general, and hybrid male sterility in particular, can have several nonmutually exclusive explanations (Coyne, 2018;Coyne & Orr, 2004), including X-linked gene misregulation, the disruption of the dosage compensation, and inactivation of the X chromosome during spermatogenesis (Good et al., 2010;Johnson & Lachance, 2012;Presgraves, 2008;Rodriguez et al., 2007). Nevertheless, the X chromosome does not always show enrichment for misregulated genes compared to the autosomes (Gomes & Civetta, 2015;Mugal et al., 2020). This picture might be compounded by the acquisition of a new X chromosome (aka neo-X), which will impact the patterns of molecular evolution of the resident genes. Depending on the time of origination of the neo-X relative to the ancestral X (anc-X) chromosome, the two chromosomes might play a different role in the divergence of expression profiles and regulatory networks as the properties of many of the genes harboured in the two chromosomes can substantially differ (Assis et al., 2012;Meisel et al., 2012b). How an anc-and a neo-X chromosome can impact expression differentiation across different reproductive tissues between recently diverged lineages and their hybrids remains largely unexplored.
The fly Drosophila willistoni has a metacentric X chromosome that resulted from the fusion between the chromosomal elements A and D of the ancestral Drosophila karyotype (Muller, 1940).
Specifically, Muller's element A corresponds to the ancestral X chromosome in the genus Drosophila whereas Muller's element D, an autosome in most Drosophila species (Table S1), is an acquired neo-sex chromosome (Spassky & Dobzhansky, 1950). D.
willistoni includes three allopatric subspecies, all of them found in the American continent: D. willistoni willistoni; D. w. winge; and D. w. quechua (Mardiros et al., 2016). The first two subspecies are not only indistinguishable at the morphological level, including fast evolving traits such as the male genitalia (Civetta & Gaudreau, 2015), but also in terms of sequence differentiation as revealed by the gene mtCOI (Mardiros et al., 2016). Although no fixed premating isolation has been found between these subspecies (Davis et al., 2020), they do exhibit unidirectional hybrid male sterility (Mardiros et al., 2016). Notably, and unlike other sterile hybrids between recently diverged Drosophila lineages that show sperm defects and reduced sperm motility (Gomes & Civetta, 2014;Kulathinal & Singh, 1998;Moehring et al., 2006), these sterile hybrid males produce typical numbers of motile sperm but cannot transfer them during copulation (Civetta & Gaudreau, 2015; Gomes & Civetta, 2014). This failed sperm transfer results from a bulging of the seminal vesicle, which hinders sperm migration towards the vas deferens, and therefore sperm incorporation into the ejaculate, which is otherwise successfully transferred to the  (Davis et al., 2020). Collectively, the features of the D. w. willistoni and D. w. winge subspecies are largely suggestive of a very early stage in their process of consolidation as distinct species, offering an excellent opportunity to detect misregulation in their hybrids while discerning the role of different regulatory elements on expression divergence under minimal biases introduced by dissimilar cell type composition. Moreover, this subspecies pair is an ideal system to analyse the role of an anc-and neo-X chromosomes during early stages of transcriptome differentiation.
Here, we have constructed a new reference-quality genome assembly using long-sequencing reads and upgraded its annotation. With these new 'omic resources, we have characterized the transcriptome of three tissues relevant for reproductive functions (testis, male accessory gland, and ovary) using RNA-seq in D. w. willistoni and D. w. winge, aiming at: (i) uncovering transcriptome differentiation at the tissue-level between parental subspecies and their sterile F1 hybrid; (ii) assessing whether the tissuespecific expression patterns of the two arms of the X chromosome have converged while testing if this chromosome shows disrupted gene dosage compensation in testis; (iii) scrutinizing the types of transgressive expression across tissues between the sterile hybrid and the parental subspecies while determining how they are distributed in the genome; (iv) deciphering the architecture of regulatory evolution behind the differentiation of the subspecies by conducting ASE analyses; and (v) extracting functional interpretable patterns from transgressively expressed genes in interspecific hybrids to determine how gene functional attributes and protein-protein interaction networks might have been impacted early in speciation. The new omic resources generated, as well as the portrait of expression differentiation between the subspecies and their hybrids, make the D. willistoni subspecies system particularly appealing for gaining key insight into the interface between transcriptome evolution and the early stages of speciation. winge Uruguay males; in this cross, hybrid males are sterile (Mardiros et al., 2016). Vials were maintained at room temperature, under 24 h light conditions, and on standard corn meal-agar-dextrose medium.

| Fly husbandry
Adult flies were collected as newly emerged within 4 h to ensure virginity. Separate sets of flies were collected: 3-5 day-old virgin flies for setting crosses or for total RNA extraction with annotation purposes; 5 day-old virgin flies for tissue-level total RNA extraction utilized in differential expression assays; and 0-4 day-old individuals for DNA extraction to be utilized for genome sequencing.

| Genomic DNA extraction and de novo genome assembly construction
Genomic DNA from females of the standard strain was extracted for Illumina paired-end (PE; 150 bp) sequencing .
Library preparation was performed using the NEXTFLEX DNA kit (insert size = 403 bp), and sequencing was done over a quarter of a lane on an Illumina HiSeq 2500 instrument. Further, genomic DNA for Single Molecule Real-Time (SMRT) sequencing from PacBio technologies was prepared from 160 2 hr-starved males per column of the same strain . Sheared genomic DNA was size selected from 20-80 kb long fragments using the BluePippin system.
The library was sequenced using PacBio's P6-C4 chemistry across 17 cells. For assembly generation, we implemented a computational pipeline that integrates Illumina PE150 and PacBio sequencing reads  (Supporting Information S1).

| Physical mapping
Fifty-two protein-coding genes were PCR-amplified (Table S2 for primer sequences), cloned, and mapped by in situ hybridization on polytene chromosomes of the standard strain (Ranz et al., 1997).
Mapping information from eight additional genes was added (Chan et al., 2015).

| RNA sequencing
Tissue dissection (testis, male accessory gland, and ovary) was performed daily during 2 h windows in ice-cold 1 × PBS and always trying to randomize (e.g., the order of dissection of different genotypes) or minimize (e.g., the time of the day the dissection was started) any possible experimental error. Total RNA extraction and quality control were performed as reported (Clifton et al., 2017;Ranz et al., 2021). Twenty-seven ribodepleted, strand-specific libraries (3 biological replicates × 3 genotypes × 3 tissues) were prepared . The resulting libraries were multiplexed and PE sequenced (100 bp) over one lane per tissue on an Illumina HiSeq 2500 instrument. An additional library was identically prepared using total RNA from whole bodies of naïve males and virgin females, and sequenced over one lane on the same instrument. All RNA and genome sequencing was performed at the UCI Genomics High-Throughput Facility.

| Repeat and gene annotation
Repeat content was predicted, combined with existing information in Diptera, and used to soft mask the genome assembly generated.
Protein-coding gene models in the release dwil_r1.05 (FB2017_04) of FlyBase (dos Santos et al., 2015) were identified using Exonerate 2.47.3 (Slater & Birney, 2005) and upgraded by leveraging the RNA-seq libraries generated. This was done using StringTie version 1.3.2d (Pertea et al., 2016) and PASA version 2.3.3 (Haas et al., 2003). For miRNAs, tRNAs, rRNAs and other types of short noncoding RNAs, we combined existing information from different databases with dedicated computational tools. Lastly, lncRNA identification required removing protein coding gene models from the StringTie merged predictions and using FEELnc (Wucher et al., 2017), demanding that lncRNA gene models should have associated >200 nt long transcripts, include ≥1 splicing junction, do not overlap with rRNA gene models, and be antisense if they overlapped protein-coding gene models (Text S2).

| Differential gene expression analysis and tissue preferential expression
RNA-sequencing reads were examined for quality control and mapped against our de novo genome assembly with STAR (Dobin et al., 2013) (Supporting Information S3). Pairwise differential expression analysis across the parental subspecies and their hybrids was performed per tissue using both edgeR (Robinson et al., 2010) and DESeq2 (Love et al., 2014). A log 2 fold-change threshold of 0.5 was further applied to increase the true positive rate (Schurch et al., 2016), and a 5% false-discovery rate (FDR) implemented (Benjamini & Hochberg, 1995). The intersect between the lists of differentially expressed genes from the two methods was used. All tools used for the differential gene expression analysis were ran on UseGalaxy (http://usega laxy.org). Further, for deeming a gene as tissue-preferentially expressed, we used several criteria: a minimum threshold expression value; a minimum fold-change in expression between the tissues compared; and a minimum expression-specificity index ≥0.9 (Supporting Information S3).

| Allele-specific expression analysis to identify cis-and trans-regulatory incompatibilities
Fixed subspecies-specific single nucleotide polymorphisms (SNPs) were identified (Supporting Information S4) using Naïve variant caller followed by the Variant annotator (Blankenberg et al., 2014), and used to estimate their relative presence in hybrids. Criteria for the identification of fixed SNPs, procedures followed to measure allele specific expression in the hybrids, and the mode of regulatory evolution for each gene were determined according to the outcome across three statistical tests implemented (Fraser, 2019;Go & Civetta, 2020;Gomes & Civetta, 2015;McManus et al., 2010).
Briefly, significant expression differences between the parental subspecies (P GUA vs. P URU ) and between their alleles in the hybrid (H PGUA vs. H PURU ) were tested using a binomial exact test. To keep consistency with the differential expression analysis, a 0.5 log 2 fold-change threshold was applied in addition to the binomial exact test. To detect significant differences between the ratio of parental expression with the ratio of each parental allele in the hybrid (P GUA /P URU vs. H PGUA /H PURU ), a Fisher's exact test was used. A 5% FDR was applied to the p-values resulting from both the binomial exact test and Fisher's exact test (Benjamini & Hochberg, 1995).

| Functional analyses
Protein-protein interactions among focal genes were predicted using STRING version 11.0 (Szklarczyk et al., 2019) under the medium and high confidence thresholds, and using all active interaction sources. Significant enrichment for gene ontology (GO) terms was determined with g:Profiler (Raudvere et al., 2019). In both cases, a 5% FDR threshold was applied (Benjamini & Hochberg, 1995).

| Statistical analyses
We used built-in statistical functions in R (R Development Core Team, 2016). In the case of chi-square tests for contingency tables, the null hypothesis was systematically tested by performing 2000 Monte Carlo simulations. In the case of the differential expression analyses performed, the lists of differentially expressed genes according to the combination of criteria indicated above resulted always to be more conservative than just the application of a joint 5% FDR across tissues (Supporting Information S5).

| De novo genome assembly and annotation
We generated Illumina PE-150 and PacBio sequencing reads for the standard D. w. willistoni strain, resulting in a ~85× and ~65× coverage, respectively. Both sequencing outputs were used to generate a de novo genome assembly (dwil_UCI; Supporting Information S1 and Figure S1), which is 5% larger than the current reference assembly (dwil_caf1), with a substantially smaller gap size, fewer scaffolds, and twice the N50 value (Table 1). This enhanced contiguity is epitomized by the longest scaffold, which corresponds to the entire chromosome 3 (~34 Mb). Further, gene-level completeness was ascertained using Benchmarking universal single-copy orthologues (BUSCO) version 2.0.1 (Simao et al., 2015), recovering 98.1% (n = 3285) of Diptera genes as single copy and in complete form (Table 1). Internal assembly validation involved examining the relative order of two gene sets, including cytological mapping information, finding no disagreement (Supporting Information S6; Table S2; Figure S2). The new genome assembly was virtually complete, substantially more contiguous than the existing one, and highly reliable in terms of internal sequence order.
We annotated repetitive and low-complexity regions in the dwil_UCI assembly. In total, 102.5 Mb were populated by repeats, with 91.5 Mb (37.01%) corresponding to interspersed repeats. Our assembly includes 8.7% more repeats (41.5% vs. 32.8%, Table S3), particularly from the long terminal repeat (LTR) class (Table S3). This is expected as a more contiguous assembly is less likely to collapse repeat sequences compared to more fragmented assemblies. Further, we upgraded the models of protein-coding genes in the annotation dwil_r1.05 by incorporating RNA-seq data primarily from the reproductive tissues (testis, accessory gland, and ovary), which are the focus of our transcriptome comparative analysis (Supporting Information S7; Tables S4 and S5). The improvements affect many gene features, including an increase in the number of transcripts per gene and number and span of UTRs and exons, resulting in for example the correction of previously truncated gene models ( Figure S3) and a more reliable quantification of expression levels in downstream analyses. The new annotation (dwil_UCI_r1.0) also increases substantially the number of lncRNAs, from 120 to 1679 (Table S4; Supporting Information S7). Overall, 28.51% of the genome (70.8 Mb) is transcribed into primary transcripts with 10.3% associated with mature transcripts, and coding sequences representing 8.88% of the genome. Lastly, we used protein-coding gene models from the two annotations to anchor the contigs of the new assembly to the Muller's chromosomal elements in D. willistoni (Supporting Information S8; Table S6).

| Transcriptome sequencing
We sequenced in triplicates the transcriptome of testis, accessory gland, and ovary of the parental subspecies D. w. willistoni and D. w. winge, and their F 1 hybrids, generating 335 million of RNA-seq reads (Table S7). We found no evidence of bias between the transcripts of the parental subspecies when mapped against the reference genome (Supporting Information S9; Table S8). Principal component analysis largely substantiated the expected grouping of samples ( Figure S4), confirming the reliability of the sequencing results.
Further, we found significant differences in the fraction of uniquely mapped reads across tissues (Kruskal-Wallis rank of sums, p < .0001), with ovary showing the highest numbers (nonparametric pair-wise Steel-Dwass test; ovary vs testis, P adj = .0012; ovary vs accessory gland, P adj = .0012; testis vs. accessory gland, P adj = .0012). This could be due to an incomplete ribosomal depletion in the accessory gland and testis samples, an unanticipated biological difference across tissues, or both (Supporting Information S10; Figure S5).
Further, neither coding nor lncRNA genes showed any bias in their patterns of differential expression towards any particular parental subspecies across tissues (3-sample test for equality of proportions; TA B L E 1 Salient features of an existing assembly for D. willistoni and the newly generated dwil_UCI. coding: χ 2 = 5.88, d.f. = 2, P adj = 0.106; lncRNA, χ 2 = 4.43, d.f. = 2, P adj = 0.109; Figure 1b). In terms of genomic distribution, and unlike other reports (Brawand et al., 2011;Llopart, 2012;Meisel et al., 2012b), we did not find evidence that the X chromosome was enriched for differentially expressed genes relative to the autosomes in any of the tissues assayed (Table S17) Note: Gua, D. w. willstoni; Uru, D. w. winge. Differentially expression was set at a 5% FDR and a 0.5fold change in expression. Direction of the differential expression between the two subspecies: > overexpression, < underexpression. a In parenthesis the number of lncRNAs, rRNAs, tRNAs, and snoRNA, respectively.
b Only genes expressed across the three types of biological samples.
c Genes that show differences in mRNA levels for at least one tissue in a given direction between the subspecies that are not observed in at least one other tissue.
TA B L E 2 Salient patterns from comparative expression analysis between D. willistoni subspecies.
F I G U R E 1 Relationship between tissue types and patterns of differential expression between the parental subspecies. (a) Percentage of nondifferentially (NDE) and differentially (DE) expressed genes between the parental subspecies. Testis harbour significantly more differentially expressed genes than accessory gland and ovary. (b) Percentage of differentially expressed genes relative to the subspecies in which they exhibit significantly higher (>) or lower (<) expression. Only coding and lncRNA genes were considered. Gua, D. w. willstoni; Uru, D. w. winge.

| Tissue-preferential expression and the demasculinization of the X chromosome
Previous studies in Drosophila documented an underrepresentation of male-biased -demasculinization-and an overrepresentation of female-biased -feminization-genes in expression on the X chromosome (Assis et al., 2012;Parisi et al., 2003;Sturgill et al., 2007). The X chromosome of D. willistoni is the result of a fusion between an anc-X (left arm or XL) and a neo-X (right arm or XR). As the pattern of masculinization and feminization of the gene content of the anc-and neosex chromosomes can differ in magnitude (Assis et al., 2012;Meisel et al., 2012b), D. willistoni offers the opportunity to test for such differences. In fact, previous analyses with this in mind were performed in this species (Meisel et al., 2012a), but involving anatomical sections (abdomen, thorax, and head), which are amalgams of different tissues.
By using our tissue-specific data and more complete gene annotation, we revisited the extent to which the anc-and neo-X chromosomes of D. willistoni differ in their degree of feminization/demasculinization.
Different types of tissue-preferentially expressed genes in the parental subspecies showed very different chromosomal distributions (χ 2 = 93.59, d.f. = 2, P adj = 3.3 × 10 −16 ; Table S18). Ovarypreferentially expressed genes are overrepresented on the anc-X but not on the neo-X (post hoc tests, P adj < .001 and P adj > .05, respectively; Table S18). Accessory gland-preferentially expressed genes were significantly underrepresented on the anc-X and overrepresented on the autosomes, particularly on the 2R, as reported in other Drosophila species (Meiklejohn & Presgraves, 2012;Mikhaylova & Nurminsky, 2011). Testis-preferentially expressed genes were not significantly depleted from the anc-X, being also present according to random expectation on the neo-X. Grouping the Muller's elements based on whether they form part of the X chromosome or the autosomal complement did not change these conclusions nor did when more restricted fold changes in expression, four-fold instead of two-fold, were applied. Lastly, we re-examined our findings by considering lncRNA and protein-coding genes separately (Table S19). These two classes of genes did not necessarily show consistent patterns of chromosome distribution. For example, testispreferentially expressed lncRNA genes were enriched in the anc-X while they were either underrepresented (2R) or found according to random expectation in the neo-X and the autosomes (2L and 3). In contrast, testis-preferentially expressed coding genes do not show any chromosomal bias. Overall, these results highlight the very different dynamics between lncRNA and coding genes in D. willistoni while uncovering notable differences with previous observations in other Drosophila species (Supporting Information S12).

| Chromosome-level mechanisms of expression in the hybrid and parental subspecies
Whole-chromosome-based regulatory mechanisms such as incomplete dosage compensation and meiotic sex chromosome inactivation, if disrupted, can contribute to sterility in interspecific hybrids (Good et al., 2010;Johnson & Lachance, 2012;Presgraves, 2008;Rodriguez et al., 2007). To evaluate this possibility, we performed two analyses (Supporting Information S13). We first ascertained whether the X:A expression ratios was lower than 1 by calculating absolute median expression levels in testis and accessory gland (X:AA) as well as in ovary (XX:AA). We found significantly lower median expression levels in the two arms of the X chromosome relative to the autosomes in testis and accessory gland. In testis, the neo-X arm was expressed 28%-32% less, depending on the genotype, than the autosomes, with the accessory gland always showing a more moderate reduction in expression (11%-14%). The anc-X shows a similar significantly reduced expression relative to the autosomes in testis (27%-30%) as well as in accessory gland (20%-26%). Conversely, in ovary, we did not detect statistically significant differences in expression among the anc-X, the neo-X, and the autosomes (Figure 2; Figure S8). Overall, we found virtually the same patterns of interchromosomal expression across three minimum expression thresholds considered ( Figures S9 and S10). This is suggestive of no overt difference in the degree of downregulation of the X-chromosome in the two tissues examined in sterile hybrid males in relation to the males of the parental subspecies.
Further, the analysis of the degree of expression equalization between the sexes, that is, that based on the comparison of ratios between the two male tissues to ovary for both the autosomes (AA:AA) and the two arms of the X (anc-X:anc-Xanc-X; neo-X:neo-Xneo-X), allowed to further assess the dose effect in the face of incomplete dosage compensation (Figures S11-S13) (Mank, 2013).
Additionally, it also revealed that the anc-X displays a significantly lower expression ratio in both the comparison testis to ovary and that of accessory gland to ovary, regardless of the minimum expression threshold considered. For the neo-X, the pattern is similar in the comparison testis to ovary while it is dependent on the minimum expression threshold considered in the case of the comparison accessory gland to ovary. These results mirror the unequal distribution of tissue-preferentially expressed genes between the two arms (Table S18), further suggesting that the neo-X chromosome is in a delayed stage of demasculinization compared to the anc-X (Supporting Information S14).

| Patterns of differential expression in hybrids across tissues
As transgressive expression can be relevant for hybrid sterility (Brill et al., 2016;Civetta, 2016;Go & Civetta, 2020;Mack & Nachman, 2017;Moehring et al., 2007), we identified those genes that were either over or underexpressed in the hybrid in relation to both parental subspecies. In addition, we also identified genes showing additive expression patterns, that is, those expressed higher in the hybrid than in one of the parental subspecies but lower relative to the other (Figure 3a; showed transgressive expression in more than one tissue (Figure 3b), underscoring that transgressive expression in the hybrids examined is effectively tissue-dependent.
Next, we analysed whether transgressively expressed genes were preferentially expressed in particular tissues as, due to their narrower expression breadth, they are thought to be less pleiotropic and under weaker selective constraints, which might result in a higher rate of evolution (Assis et al., 2012;Khaitovich et al., 2005;Meisel, 2011). Under this assumption, we expect that tissuepreferentially expressed genes in the parental subspecies would be more likely to show transgressive expression in the hybrids. We did not find a prevalence of tissue-preferentially expressed genes among those transgressively expressed, with testis showing even a lower proportion than ovary (Supporting Information S15 and Table S20). Further, to discard any transgressive expression in testis as the result of differential cell-type composition in the hybrid (Good et al., 2010;Mugal et al., 2020), we examined the expression levels of 10 known cell types marker genes (Witt et al., 2019;Witt, Shao, et al., 2021) across the hybrid and parental subspecies. None of the genes shows statistically significant differential expression between the parental subspecies ( Figure S14). More importantly, only one F I G U R E 2 Absolute expression levels of the XL, XR, and the autosomes in three tissues and three genotypes of D. willistoni. The violin plots and nested box plots show the distribution of the log 2 (expression values) for genes expressed >1 TPM (transcripts-per-million). Autosomes are indicated in blue while the ancestral and neo-X chromosomes are in orange (XL) and red (XR), respectively. Medians are indicated on top of each plot. Statistical significance of differences between median values was determined using Kruskal-Wallis tests, followed by pairwise comparisons using Wilcoxon rank sum tests and correcting for multiple tests, which was followed by an additional correction that considered the number of tests performed (Benjamini & Hochberg, 1995). Solid lines indicate particular chromosomal contrasts entailing a statistically significant difference in expression: P adj < .05 (*); P adj < .01 (**); and P adj < .001 (***). The number of genes considered is shown on the top left of each plot. Gua, D. w. willstoni; Uru, D. w. winge; Hyb, hybrid resulting from a Gua female and Uru male. ** ***

| Genomic distribution of hybrid differential expression
The disproportionate role of the X chromosome to hybrid disfunction through gene misregulation has been documented in some but not all cases (Gomes & Civetta, 2015;Llopart, 2012;Masly & Presgraves, 2007;Mugal et al., 2020). We examined whether differentially expressed genes in the hybrids relative to the parental subspecies were equally abundant on the X chromosome and autosomes, and whether there was any difference among tissues. For 493 genes with reliably determined chromosomal location, we found an uneven distribution of the additive and transgressive patterns in testis but not in ovary or accessory gland, in this last case possibly because of limited statistical power ( Table 3). In testis, the two arms of the X chromosome seemed enriched for transgressively expressed genes, while the autosomes were enriched for additively expressed genes (χ 2 = 32.68, P adj = 1.5 × 10 −3 ). The two types of transgressively expressed genes were not found to be differentially associated at the chromosomal level (χ 2 = 0.778, p = .752).
We also sought for evidence of clustering of transgressively expressed genes at a finer scale, which could denote misregulation Note: P adj values are provided in parenthesis next to the trend symbol; ↑, enrichment; ↓, depletion; ns, no deviation in relation to the random expectation. P adj values correspond to the implementation of the Benjamini-Hochberg correction (Benjamini & Hochberg, 1995). a As the chi-square test cannot be computed with zero marginals, a G-test of independence was performed instead.
of a regulatory element of regional effect or differential chromatinbased regulation (Degner et al., 2012;McVicker et al., 2013). We used three arbitrary cutoff values (5, 10, and 15 kb) to determine the share of genes dubbed as clustered, as well as the number of clusters. We found abundant clustering at the three cutoff values (Table S21). Notably, 12 clusters remained detected across cutoff values, seven of them being relevant for ovary, and another five for testis. When examined these clusters for the presence of gene models belonging to the same orthogroup, we found that two of the clusters relevant for ovary included paralogues.

| Identification of genome-wide regulatory mismatches in interspecific hybrids
Divergent regulatory elements in a hybrid genome can cause differential expression associated with the different effects of the parental alleles (Coolon et al., 2015;Landry et al., 2005;McManus et al., 2010;Wittkopp et al., 2004). To determine the extent of such regulatory mismatches, we used putatively fixed SNPs between the parental subspecies to identify allele-specific gene expression in the hybrids. Upon identifying usable SNPs for the three tissues assayed (Table S22), we inferred modes of regulatory divergence (Materials and methods; Table 4). As this requires the presence of both parental alleles in the hybrid, the analysis could only be performed on autosomal genes for testis and accessory gland. For ovary, the analysis was done separately for the autosomal and X-linked genes. Altogether, we could investigate the mode of regulatory divergence for 1626 genes in testis, for 1821 in accessory gland, and for 6384 in ovary.
Regardless of the tissue, 81%-88% of the analysable genes showed no evidence of regulatory divergence (Table 4). Of the remaining genes, regulatory divergence in the gonads (i.e., testis and ovary) appear to be more significantly impacted by cis-than by trans-changes, with no significant differences being observed in accessory gland (Table S23). Notably, 89.8% of the genes analysable in at least two tissues showed consistency in their mode of regulatory evolution, predominantly in association with the conserved category. Further, we detected very little evidence of compensatory evolution as well as very limited presence of artefactual significant negative correlations between cis and trans effects (Supporting Information S16; Table S24; Fraser, 2019;Zhang & Emerson, 2019).
We then analysed the relationship between the mode of regulatory evolution and transgressive expression in the hybrids, which was possible for testis and ovary but not for accessory gland due to the limited number of genes showing transgressive and additive expression patterns (five and four genes, respectively). For both tissues, we found evidence of nonrandom associations between particular modes of regulatory evolution and patterns of differential a Different modes of regulatory divergence as inferred from the allele specific expression analysis, which entails the comparison of expression levels between the parental subspecies (P GUA vs. P URU ), the comparison of their alleles in the hybrid progeny (H PGUA vs. H PURU ), and the comparison of the expression ratios between the parental alleles in parental individuals and in the hybrids (P GUA /P URU vs. H PGUA /H PURU ). Six modes of regulatory divergence are considered. Conserved refers to genes for which there is no evidence of regulatory evolution in either cis or trans. No significant differential expression is observed between parental subspecies (P GUA = P URU ), between their alleles in the hybrids (H PGUA = H PURU ), and between the ratio of parental alleles expression and the ratio between the parental alleles in the hybrids (P GUA /P URU = H PGUA /H PURU ). Cis only refers to genes whose regulatory evolution is driven by genetic differences in their cis regulatory elements. Significant differential expression is detected between the parental subspecies and between their alleles in the hybrid but no significant differences are detected between the ratio of parental alleles expression and the ratio of the parental alleles in the hybrid. Trans only refers to genes whose expression evolution is driven by changes affecting the trans-acting factors that regulate them. Significant differential expression is detected between the parental subspecies but not between the parental alleles in the hybrid. Significant differences between the ratio of parental allele expression and that of the parental alleles in the hybrid are detected.
Cis & trans refers to genes that are evolving at the regulatory level through changes at both levels. Significant differences in expression are detected between parental subspecies, in the expression of the parental alleles in the hybrid, and between the ratio of parental allele expression and the ratio of the parental alleles in the hybrid. Compensatory refers to genes that show no significant expression differences between the parental subspecies despite having evidence of both cis and trans divergence. Regulatory divergence in both cis-and trans-perfectly compensate for each other, resulting in no significant differential expression between the parental subspecies. Significant differences in expression between alleles in the hybrid and between the ratio of parental expression and the ratio of alleles in the hybrid are nevertheless detected. Ambiguous refers to genes for which the expression patterns observed in the parental subspecies and the hybrid do not fall into any of the above patterns. Since males are hemizygous for the X chromosome, only autosomal genes were available for analysis in testis and accessory gland samples. For ovary, both autosomal and X-linked are analysable (Table S22). Only comparable results (i.e., autosomal genes) across tissues are shown.
expression. Bootstrapping showed that in ovary, but not testis, genes experiencing cis-regulatory divergence exhibit the most significant enrichment among those showing additive expression (Figure 4). In contrast, genes displaying transgressive expression, both in ovary and testis, were primarily characterized by conserved regulatory evolution ( Figure 4).
Lastly, we tested the expectation of a comparatively larger effect on expression of cis relative to trans changes due to the presumed more limited detrimental pleiotropic effects of the former (Coolon et al., 2014;Hill et al., 2021;Kopania et al., 2022). For testis, the size effect of the expression difference induced by cis changes is significantly larger than that of trans changes, while for accessory gland and ovaries we found no significant difference (Wilcoxon rank sum test: testis, P adj = 0.0429; accessory gland, P adj = 0.536; ovary, P adj = 0.536; Table S25). This lack of consistency across tissues might reflect differences in the strength of purifying selection operating on them. Collectively, and despite the prevalence of the conserved mode of regulatory evolution across the tissues of these two D.
willistoni subspecies, subtle differences in the role and magnitude effects of cis-and trans-regulatory changes, as well as different associations to additive and transgressive expression across tissues, highlight the importance of the tissue context to properly uncover patterns of regulatory evolution.

| Functional patterns of enrichments among differentially expressed genes in interspecific hybrids
Cis-and trans-regulatory incompatibilities in a hybrid background could lead to a cascade of differential expression (Go & Civetta, 2020), being preferentially associated with particular functional rubrics or part of a shared network of interactions. We examined these two features (Materials and methods) assuming functional conservation between D. willistoni and D. melanogaster for the genes examined.
For accessory gland, neither significant protein-protein interaction (PPI) network nor overrepresentation of biological process were detected. In ovary, 55 out of the 228 differentially expressed genes in the hybrids were predicted to have more protein-protein interactions than randomly expected (PPI enrichment P adj = 4.34 × 10 −13 ).
In addition, this set of genes were disproportionally enriched for the GO biological process term chorion-containing eggshell formation (P adj = 0.0307). In testis, we also found evidence for a significant PPI network involving 64 out of the 263 differentially expressed genes (PPI enrichment P adj = 8.11 × 10 −9 ), and up to 12 functionally enriched biological process terms, including cell adhesion (P adj = 0.0307), which has been also found enriched in the transcriptome analysis of sterile male hybrids between two other very closely related Drosophila subspecies (Go et al., 2019;Gomes & Civetta, 2015) ( Table 5; Figures S15-S16; Supporting Information S17). Restricting these analyses to the 122 and 177 transgressively expressed genes in ovary and testis, respectively, does not essentially alter the patterns of enrichment found for biological processes, and many shared interaction networks remained ( Figure 5; Tables 5 and S26). A qualitative difference between these most restricted interaction networks in ovary and testis is the prevalence of transgressive overexpression in the latter. Notably, in the case of transgressively expressed genes that are part of PPI networks, we found that seven of them were part of three physical gene clusters.

| DISCUSS ION
We have comprehensively compared the transcriptome of two D.
willistoni subspecies, D. w. willistoni and D. w. winge, and the hybrids resulting from the cross between them that yields sterile males F I G U R E 4 Nonrandom association between regulatory divergence between subspecies and differential expression in the hybrid in relation to those subspecies. Heat maps for 53 and 71 genes differentially expressed in testis and ovary between the hybrid and the parental subspecies, and in relation to different types of regulatory evolution, are shown. The category of ambiguous expression divergence is excluded, and those corresponding to "compensatory" and "cis and trans" are combined and labelled as "cis + trans". Bootstrapping (n = 1 × 10 5 ) was used to assess departures from the random expectation: p < .05, one arrowhead; p < .01, two arrowheads; and p < .001, three arrowheads. The directionality of the arrow denotes whether the significant association represents enrichment or depletion. (Davis et al., 2020;Mardiros et al., 2016). This unidirectional hybrid male sterility plus a negligible cell-type tissue composition differentiation in the sterile hybrid make this subspecies pair well poised to address how regulation in gene expression evolves during the early stages of speciation across tissues, resulting in hybrid dysfunction.
Further, we have generated a more contiguous genome assembly, which should facilitate the study of the rampant structural variation in the D. willistoni species group (Rohde & Valente, 2012), and an TA B L E 5 Patterns of enrichment of gene ontology terms among genes differentially expressed in hybrids in relation to the parental subspecies. ; redundant terms were excluded using REVIGO (Supek et al., 2011). P adj values correspond to the implementation of the Benjamini-Hochberg correction to the number of terms and tissues considered (Benjamini & Hochberg, 1995). a It includes both transgressive and additive patterns of differential expression in the hybrid relative to the parental subspecies.
F I G U R E 5 STRING protein-protein interactions (PPI) networks for transgressively expressed genes in the testis and ovary of hybrids relative to the parental subspecies. Lines between nodes denote interacting proteins. Only high-confidence interactions are considered. PPI enrichment: testis, P adj = 4.44 × 10 −08 ; ovary, P adj <3 × 10 −16 . A colour code identifies different types of differential expression: red, overexpression; green, underexpression. Purple stars, cell adhesion related genes. Blue stars, chorion-containing eggshell formation related genes.
upgraded annotation with more reliable gene models and a substantially higher number of lncRNA genes, which are relevant both in transcriptional regulation and phenotypic change (Wen et al., 2016;Zhu et al., 2021).

| Expression profile differences across tissues in D. willistoni
Not many studies have focused on the early lineage divergence, at the onset of hybrid breakdown, comparing gene expression across tissues (Mugal et al., 2020). Prior comparative expression analyses among male-, female-, and nonbiased genes have shown that the former category evolves the fastest in Drosophila (Llopart, 2012;Ranz et al., 2003). Here, we also found that testis exhibit not only a higher fraction of differentially expressed genes between subspecies than ovary but also in relation to the male accessory gland, which although somatic, it is also part of the reproductive system.
Roughly, one fourth of all the genes consistently expressed across tissues are differentially expressed between D. w. willistoni and D.
w. winge in at least one of the tissues assayed. Notably, the degree of consistency in relation to the identity of the genes that are differentially expressed between the subspecies was minimal, stressing the need of tissue-level characterization even during early lineage differentiation.

| The X chromosome arms show very distinctive patterns of gene content differentiation
Early results pointed to the demasculinization of the gene content of the X chromosome relative to the autosomes (Assis et al., 2012;Parisi et al., 2003;Sturgill et al., 2007). A key explanation is that the X chromosome spends more time in females than in males, which should result in an overrepresentation of genes beneficial for females at the expenses of males (Rice, 1984;Vicoso & Charlesworth, 2006). Our results reveal that the X chromosome of D. willistoni is not depleted of testis-preferentially expressed genes. In fact, when expression breadth is considered, we find no evidence of depletion of testis-specific genes on the X chromosome as previously shown in D. melanogaster (Meiklejohn & Presgraves, 2012;Meisel et al., 2012a). Further, when considering lncRNA and protein-coding genes separately, the former are overrepresented on the anc-X while the second follow the random expectation. This differs from observations on the X of D. pseudoobscura, species in which lncRNAs follow the random expectation (Nyberg & Machado, 2016). The reason for this interspecific difference is not apparent currently.
The anc-and neo-X chromosomes of D. willistoni exhibit additional differences. The anc-X shows depletion of accessory gland-preferentially expressed genes (demasculinization) and an overrepresentation of ovary-preferentially expressed genes (feminization). These findings agree with previous observations on the X of D. melanogaster (Meiklejohn & Presgraves, 2012;Meisel et al., 2012a). The neo-X chromosome does not deviate from the random expectation for these two gene classes. Collectively, the anc-and neo-X chromosomes of D. willistoni are characterized by marked differences in gene content composition, consistent with a neo-X chromosome that has not yet converged functionally with the anc-X chromosome, similarly to D. pseudoobscura (Assis et al., 2012).

| No evidence of disrupted dosage compensation in male sterile hybrids
In several Drosophila species, lower X-to-autosome expression ratio in testis and accessory gland has been reported (Assis et al., 2012;Meiklejohn et al., 2011;Meisel et al., 2012a;Parisi et al., 2003) while a higher X-to-autosome expression ratio has been found in D. melanogaster ovary (Assis et al., 2012;Parisi et al., 2003) although not in D. pseudoobscura (Assis et al., 2012). The patterns detected here for D. willistoni mirrored those observed in D. pseudoobscura for all the tissues and genotypes assayed.
In the case of testis, a lower X-to-autosome expression ratio could result from gene silencing on the X chromosome as spermatogenesis progresses (Mahadevaraju et al., 2021), lack of dosage compensation (Meiklejohn et al., 2011), or demasculinization of its gene content due to intralocus sexual antagonism (Charlesworth et al., 1987;Rice, 1984;Vicoso & Charlesworth, 2006). This last mechanism though does not seem to be that relevant in D. willistoni as we do not find a depletion of testis-preferentially expressed genes on the X. Assuming similar mechanisms between D. willistoni and D.
melanogaster, canonical dosage compensation could be operating in spermatogonia and whole X-chromosome silencing, or incomplete noncanonical dosage compensation, could be operating in spermatocytes (Gupta et al., 2006;Mahadevaraju et al., 2021). Notably, the reduced gene expression on the X chromosome in D. willistoni testis is not as pronounced as in D. melanogaster (Meiklejohn & Presgraves, 2012). While the reasons for this interspecific discrepancy remain unclear, both in hybrids and parental subspecies, our results support an overall downregulation of the X chromosome during spermatogenesis -although not to the extent of D. melanogasterand a reduced X chromosome compensation in male accessory gland similarly to D. melanogaster (Meisel et al., 2012a). Therefore, dosage compensation mechanism is not disrupted in reproductive tissues of D. willistoni sterile hybrid males.

| Distinct patterns of hybrid misexpression across tissues and gene ontologies
Transgressive expression varied substantially across tissues in D.
willistoni hybrids. Testis exhibited the largest proportion of transgressively expressed genes. Interestingly, testis showed the lowest proportion of tissue-preferentially expressed genes among those with transgressive expression. While spermatogenesis genes are commonly misregulated in sterile hybrids Li et al., 2016;Moehring et al., 2007), we did not find them overrepresented among misregulated genes in hybrids. This result agrees well with the high numbers of motile sperm produced by the sterile hybrid males (Gomes & Civetta, 2014).
The pattern of gene ontology enrichment among the transgressively expressed testis genes did not reveal an overt relationship with testis function. Instead, we found several GO terms related to different aspects of cell adhesion, which is important for the proper assembly of the muscle sheaths in the D. melanogaster male reproductive system (Kuckwa et al., 2016;Susic-Jung et al., 2012). The low proportion of misregulated accessory gland genes in hybrid males aligns with the absence of differences in seminal fluid composition of the sterile males, as these are able to induce the typical changes in uterus morphology (Davis et al., 2020). Lastly, in ovary, chorion-containing eggshell formation genes were over-represented among those showing transgressive expression, possibly reflecting a fast regulatory divergence between subspecies, similarly to their evolutionary rate at the protein sequence level (Jagadeeshan & Singh, 2007).

| Evidence of a faster-X effect on hybrid differential expression
Prior expression studies in hybrids of the D. melanogaster species group, as well in hybrids between D. pseudoobscura subspecies, found fewer than expected X-linked chromosome genes misexpressed (Gomes & Civetta, 2015;Moehring et al., 2007). In contrast, sterile male hybrids in mice showed predominant overexpression of X-linked spermatogenesis genes (Good et al., 2010) and sterile male hybrids between closely related species of flycatchers had limited over and under-expression of Z-linked testis expressed genes (Mugal et al., 2020). We found a greater proportion of X-linked genes with transgressive expression in hybrids for testis than for other tissues.
It is unclear why X-linked testis-expressed genes show enrichment for transgressive expression in hybrids, but previous work has shown faster X evolution effects for male-biased and tissuerestricted genes in expression in adults (Llopart, 2012;Meisel et al., 2012b). Differences in the number of regulatory elements affecting gene expression in testis could explain differences in rate of evolution relative to other tissues. In fact, gene expression in ovary is more heavily dependent on trans-regulatory factors than in testis, which relies additionally on open promoter chromatin state (Witt, Svetec, et al., 2021). Faster rate of change in X-linked regulatory elements that preferentially regulate X-linked genes could cause a higher expression divergence of these genes relative to those on the autosomes (Coolon et al., 2015). However, we do not find the X chromosome enriched relative to the autosomes for differentially expressed genes between the subspecies, which argues against an elevated presence of cis-regulatory changes within the X chromosome. Divergence of autosomal trans-acting factors preferentially regulating the expression of X-linked genes could therefore explain our observations. Interestingly, the contribution of trans-regulatory changes to faster-X expression changes is more common within than between species in the D. melanogaster species group (Coolon et al., 2015), which suggests that the two D. willistoni subspecies analysed here mimic more closely intraspecific patterns of regulatory differentiation. Lineage-specific effects across phylogenies, developmental or age stage-specific sample effects, and whole individuals vs. tissue sampling can explain discrepancies in results among different sets of species.
Notably, we find evidence of clusters of misexpressed genes in hybrids for testis and ovary, which is suggestive of an incipient subspecies differentiation at the chromatin structure level, impacting sets of nearby genes. Although tandem duplicates are more frequently colocalized on the same topologically chromatin domain (Ibn-Salem et al., 2017;Makova & Li, 2003), we found very little evidence for their presence in the clusters identified.

| Cis-regulatory changes of variable size effects across tissues prevail during incipient speciation stages
Previous analysis on the tempo of accumulation of cis-and transacting regulatory changes during lineage divergence supported the prevalence of the former between species (Coolon et al., 2014;Gomes & Civetta, 2015;Metzger et al., 2016;Tirosh et al., 2009).
Our results in testis and ovary agree well with this observation; for accessory glands, we did not find significant differences in the occurrence of cis and trans changes. Previous work has suggested that expression changes involving diffusible trans-regulatory factors encoding genes could be buffered by selection due to possible pleiotropic effects (Lemos et al., 2008;Prud'homme et al., 2007). Further, genes expressed in testis and diverging in cis showed a larger average size effect on expression differentiation between D. willistoni subspecies than genes diverging in trans, which is compatible with more intense selective constrains acting on trans-regulatory changes. This difference was not observed in ovary and accessory gland. Overall, the preponderance of cis relative to trans changes for testis and ovary but not for accessory gland, and the different size effect that these types of changes have on expression differences between the subspecies, highlight how regulatory evolution differs across reproductive tissues. The lack of constrain observed in the accumulation of trans-relative to cis-regulatory changes in the accessory gland agrees well the majority of seminal fluid protein encoding genes evolving interspecifically under relaxed selection (Patlar et al., 2021).
A common outcome in several ASE studies has been the prevalence of compensatory changes, that is, mutations in cis and trans that differ in the direction of their effect on gene expression, thus resulting in no difference between the populations or species compared (Fear et al., 2016;Goncalves et al., 2012;Mack et al., 2016;Tirosh et al., 2009). This type of changes is expected under stabilizing selection, which preserves expression levels despite regulatory differentiation (Signor & Nuzhdin, 2018). Nevertheless, the role of compensatory evolution in ASE assays might have been overstated (Fraser, 2019;Zhang & Emerson, 2019) (but see Signor & Nuzhdin, 2019). Between the D. willistoni subspecies, we found scarce evidence of compensatory changes, which is compatible with the limited divergence time between the subspecies as well as with the high proportion of conserved genes found at the regulatory level.
This differs from other intraspecific analyses in mice (Goncalves et al., 2012), D. melanogaster (Fear et al., 2016), and during cotton domestication (Bao et al., 2019). Differences in the rate of deleterious mutations or in the efficacy of selection to eliminate them between D. willistoni and those other species could contribute to the observed difference (Andolfatto et al., 2011).

| Limitations and future directions
Our work highlights the need of investigating tissue-pattern expression when characterizing the transcriptome of interspecific hybrids while controlling for cell-type composition differences. Evidence for such differences should prompt the adoption of methodological approaches such as single-cell RNA-seq .
Further, our study is not exempt of limitations. Both the expression changes and nucleotide differences used in the ASE analysis are not necessarily fixed, as we included one line per subspecies.
Additionally, we were prone to overestimate the true divergence.
Conversely, the combination of working with early diverged lineages and only one line per subspecies can also result in genes that harbour no informative difference, limiting the power of some analyses.
Lastly, the genome sequencing of these two subspecies will allow to directly analyse nucleotide changes impacting regulatory elements as well as to identify structural changes responsible for expression differences unrelated to regulatory evolution. Such 'omics resources, plus the ones generated here, will enhance the study of speciation in the D. willistoni subgroup and across the Drosophila genus.

ACK N O WLE D G E M ENTS
We thank Richard Meisel and John Parsch for comments on the manuscript. We also thank the University of California Irvine High Performance Computing cluster for facilitating our computational analyses. This work was supported by two National Science and an NSERC Julie Payette PGS award to S.G.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All raw sequencing data have been deposited as part of the NCBI BioProject PRJNA670571. The genomic assembly UCI-dwil_1. Information S1-S4.

B EN EFIT-S H A R I N G S TATEM ENT
Benefits from this research accrue from the sharing of our data and results on public databases as described above.