
==== Front
PLoS Genet
PLoS Genet
plos
PLOS Genetics
1553-7390
1553-7404
Public Library of Science San Francisco, CA USA

39178312
10.1371/journal.pgen.1011257
PGENETICS-D-24-00414
Research Article
Biology and Life Sciences
Genetics
Gene Expression
Biology and Life Sciences
Genetics
Gene Expression
Gene Regulation
Biology and Life Sciences
Genetics
Biology and Life Sciences
Evolutionary Biology
Population Genetics
Genetic Polymorphism
Biology and Life Sciences
Genetics
Population Genetics
Genetic Polymorphism
Biology and Life Sciences
Population Biology
Population Genetics
Genetic Polymorphism
Biology and Life Sciences
Genetics
Genetic Dominance
Research and Analysis Methods
Animal Studies
Experimental Organism Systems
Model Organisms
Drosophila Melanogaster
Research and Analysis Methods
Model Organisms
Drosophila Melanogaster
Research and Analysis Methods
Animal Studies
Experimental Organism Systems
Animal Models
Drosophila Melanogaster
Biology and Life Sciences
Zoology
Entomology
Insects
Drosophila
Drosophila Melanogaster
Biology and Life Sciences
Organisms
Eukaryota
Animals
Invertebrates
Arthropoda
Insects
Drosophila
Drosophila Melanogaster
Biology and Life Sciences
Zoology
Animals
Invertebrates
Arthropoda
Insects
Drosophila
Drosophila Melanogaster
Biology and Life Sciences
Organisms
Bacteria
Gut Bacteria
Biology and Life Sciences
Genetics
Microbial Genetics
Pervasive tissue-, genetic background-, and allele-specific gene expression effects in Drosophila melanogaster
Tissue- and allele-specific gene expression effects
https://orcid.org/0000-0002-1322-1000
Glaser-Schmitt Amanda Conceptualization Formal analysis Investigation Methodology Visualization Writing – original draft Writing – review & editing 1 *
Lemoine Marion Formal analysis Investigation Visualization Writing – review & editing 2
Kaltenpoth Martin Funding acquisition Methodology Writing – review & editing 2
https://orcid.org/0000-0001-9068-5549
Parsch John Conceptualization Funding acquisition Methodology Writing – review & editing 1
1 Division of Evolutionary Biology, Faculty of Biology, Ludwig-Maximilians-Universität München, Munich, Germany
2 Department of Insect Symbiosis, Max-Planck-Institute for Chemical Ecology, Jena, Germany
Mackay Trudy F. C. Editor
Clemson University, UNITED STATES OF AMERICA
The authors have declared that no competing interests exist.

* E-mail: glaser@bio.lmu.de (AGS)
23 8 2024
8 2024
20 8 e101125714 4 2024
30 7 2024
© 2024 Glaser-Schmitt et al
2024
Glaser-Schmitt et al
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

The pervasiveness of gene expression variation and its contribution to phenotypic variation and evolution is well known. This gene expression variation is context dependent, with differences in regulatory architecture often associated with intrinsic and environmental factors, and is modulated by regulatory elements that can act in cis (linked) or in trans (unlinked) relative to the genes they affect. So far, little is known about how this genetic variation affects the evolution of regulatory architecture among closely related tissues during population divergence. To address this question, we analyzed gene expression in the midgut, hindgut, and Malpighian tubule as well as microbiome composition in the two gut tissues in four Drosophila melanogaster strains and their F1 hybrids from two divergent populations: one from the derived, European range and one from the ancestral, African range. In both the transcriptome and microbiome data, we detected extensive tissue- and genetic background-specific effects, including effects of genetic background on overall tissue specificity. Tissue-specific effects were typically stronger than genetic background-specific effects, although the two gut tissues were not more similar to each other than to the Malpighian tubules. An examination of allele specific expression revealed that, while both cis and trans effects were more tissue-specific in genes expressed differentially between populations than genes with conserved expression, trans effects were more tissue-specific than cis effects. Despite there being highly variable regulatory architecture, this observation was robust across tissues and genetic backgrounds, suggesting that the expression of trans variation can be spatially fine-tuned as well as or better than cis variation during population divergence and yielding new insights into cis and trans regulatory evolution.

Author summary

Genetic variants regulating gene expression can act in cis (linked) or in trans (unlinked) relative to the genes they affect and are thought to be important during adaptation because they can spatially and temporally fine-tune gene expression. In this study, we used the fruit fly Drosophila melanogaster to compare gene expression between inbred parental strains and their offspring in order to characterize the basis of gene expression regulation and inheritance. We examined gene expression in three tissues (midgut, hindgut, and Malpighian tubule) and four genetic backgrounds stemming from Europe and the ancestral range in Africa. Additionally, we characterized the bacterial community composition in the two gut tissues. We detected extensive tissue- and genetic background-specific effects on gene expression and bacterial community composition, although tissue-specific effects were typically stronger than genetic background effects. Genes with cis and trans regulatory effects were more tissue-specific than genes with conserved expression, while those with trans effects were more tissue-specific than those with cis effects. These results suggest that the expression of trans variation can be spatially fine-tuned as well as (or better than) cis variation as populations diverge from one another. Our study yields novel insight into the genetic basis of gene regulatory evolution.

http://dx.doi.org/10.13039/501100001659 Deutsche Forschungsgemeinschaft 274388701 https://orcid.org/0000-0001-9068-5549
Parsch John http://dx.doi.org/10.13039/501100001659 Deutsche Forschungsgemeinschaft KA2846/5-1, project number 347368302 Kaltenpoth Martin This work was supported by a Deutsche Forschungsgemeinschaft (DFG, www.dfg.de) grant (number 274388701) to JP as part of the priority program “SPP 1819: Rapid evolutionary adaptation”, and a DFG grant (KA2846/5-1, project number 347368302) to MK as part of FOR2682 “Seasonal temperature acclimation in Drosophila”. The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. PLOS Publication Stagevor-update-to-uncorrected-proof
Publication Update2024-09-05
Data AvailabilityThe gut RNA-seq and microbiome data that support the findings of this study are publicly available from National Center for Biotechnology Information (NCBI) under the BioProject accession number PRJNA1094401. Gut non-ASE gene count and expression data are available from NCBI Gene Expression Omnibus (GEO) under the series accession number GSE263264. Malpighian tubule RNA-seq data are available under the NCBI GEO accession number GSE103645. All other relevant data are within the manuscript and its Supporting Information files.
Data Availability

The gut RNA-seq and microbiome data that support the findings of this study are publicly available from National Center for Biotechnology Information (NCBI) under the BioProject accession number PRJNA1094401. Gut non-ASE gene count and expression data are available from NCBI Gene Expression Omnibus (GEO) under the series accession number GSE263264. Malpighian tubule RNA-seq data are available under the NCBI GEO accession number GSE103645. All other relevant data are within the manuscript and its Supporting Information files.
==== Body
pmcIntroduction

Gene expression variation is extensive at all organismal levels, including among tissues [1–2], cells [3–4], or alleles [5–6] of the same individual, and underlies much of the phenotypic variation that we see among individuals, populations, and species [7–9]. A long-standing challenge in evolutionary genetics has been to identify and characterize this variation. Indeed, elucidating the scope and architecture of gene expression variation as well as the mechanisms that shape it is an integral part of better understanding complex phenotypic traits [10–12], such as body size or disease susceptibility, and their evolution.

At the DNA sequence level, genetically heritable variants can modulate expression in two general ways: cis-regulatory variants, such as those within enhancers or promoters, affect the expression of linked, nearby genes, while trans-regulatory variants, such as those affecting transcription factors or regulatory RNAs, affect the expression of unlinked genes that can be located anywhere in the genome (reviewed in [13–14]). One way to interrogate the relative contribution of these types of regulatory variants to gene expression variation in species such as Drosophila, where inbred, relatively isogenic strains are available, is to compare gene expression of two parental strains or species as well as expression of their alleles in F1 hybrids [15]. Due to linkage with the allele they regulate, cis-regulatory variants affect only one of the two F1 hybrid alleles, leading to allele-specific expression (ASE), while trans-regulatory variants equally affect both alleles in the hybrid and do not lead to ASE. While cis-regulatory variation is thought to accumulate and become more predominant over larger evolutionary distances, i.e. between species [16–18], trans-regulatory variation tends to be more common among individuals within a species [5–6,19]. However, deviations from this pattern of regulatory variation have been documented in Drosophila [20–23] as well as other species [24–25], which underscores that there remains much to learn about the evolution of gene expression regulation, especially over short evolutionary distances.

An advantage of utilizing ASE to investigate the regulation of gene expression is that both the genetic basis of expression variation (e.g. cis versus trans) and the mode of expression inheritance (e.g., dominance versus additivity) can be assessed. Indeed, previous studies of ASE in Drosophila utilizing expression in F1 hybrids have found that environment [6,26], sex [27–28], genetic background [19,21,28], and body part or tissue [21–22,28] can affect regulatory architecture. However, previous studies have largely focused on single populations, long term lab strains, or comparatively closely related populations [5–6,20–22,26–28] (for an exception see [19]). Moreover, previous studies measured expression in whole animals, body parts (e.g. heads), single tissues, and/or highly functionally diverged tissues (i.e. testes versus ovaries or heads); thus, little is known about how regulatory architecture and inheritance vary among individual tissues that are spatially and/or functionally proximate. To investigate the effect of natural genetic variation from divergent populations on regulatory architecture in multiple functionally related, interconnected tissues, we analyzed messenger RNA-sequencing (mRNA-seq) data of midgut, hindgut, and Malpighian tubule tissues in four D. melanogaster strains and their F1 hybrids. Two of the strains were from a population in Umeå, Sweden [29], representing the northern edge of the species’ derived distribution, while the other two strains were from a population in Siavonga, Zambia, representing the species’ inferred ancestral range [30]. Since their divergence from ancestral populations ~12,000 years ago [31], derived D. melanogaster populations have had to adapt to new habitats, and previous studies have found evidence that at least some of the expression divergence detected between derived and ancestral African populations is adaptive [32–36].

The midgut, hindgut, and Malpighian tubules, which are analogous to the mammalian small and large intestines and kidneys, respectively, physically connect to and interact with one another at the midgut-hindgut junction and are part of the D. melanogaster digestive tract (midgut and hindgut, together with the foregut) and excretory system (hindgut and Malpighian tubules). Both systems play important roles in the regulation of homeostasis as well as the immune response [37–38] and the investigated tissues are known to engage in interorgan communication with each other, as well as with other tissues [37–39]. The excretory system is involved in waste excretion as well as ionic- and osmoregulation [38], while the digestive tract is an important modulator of food intake, nutrient absorption, energy homeostasis, and insulin secretion that can shape physiology and behavior through its interaction with the microbiome [37,40]. To investigate the effect of natural genetic variation from divergent populations on digestive tract microbiome composition, we further performed microbiome sequencing on the same gut samples for which we performed mRNA-seq. In both the mRNA-seq and microbiome data, we found extensive tissue- and genetic background-specific effects. From the ASE data, we found that genes with both cis and trans effects were more tissue-specific than genes with no differential expression regulation, although trans effects were more tissue-specific than cis effects. Despite the context specificity that we detected for regulatory architecture across tissues and genetic backgrounds, the increased specificity of trans effects was consistent, suggesting that trans-regulatory variation can be spatially fine-tuned as well as or, potentially, better than cis-regulatory variation.

Results

We performed mRNA-seq in the midgut and hindgut of two isofemale D. melanogaster strains from the northern limit of the derived species range in Sweden (SU26 and SU58) and two strains from the ancestral species range in Zambia (ZI418 and ZI197) as well as F1 hybrids between the Swedish and Zambian strains (SU26xZI418, SU26xZI197, SU58xZI418, and SU58xZI197). We additionally reanalyzed previously published mRNA-seq data from the Malpighian tubule [19] in a subset of these genotypes (SU26, SU58, ZI418, SU26xZI418, and SU58xZI418). We detected 7,675–8,209 genes as expressed in the individual tissues, with 6,894 genes that could be analyzed in all genotypes in all tissues. We focus on the genes that could be analyzed in all examined genotypes and tissues unless otherwise indicated. When considering gene expression variation across all samples, biological replicates clustered strongly by tissue type (Fig 1A). Within tissues, replicates mostly clustered by genotype, although in the hindgut there was some overlap between SU58, ZI418 and their F1 hybrid, as well as SU26 and one of its F1 hybrids (Fig 1).

10.1371/journal.pgen.1011257.g001 Fig 1 Principal component analysis of gene expression profiles in A) all examined tissues and the B) hindgut, C) midgut, and D) Malpighian tubule using all genes that could be analyzed in all or each tissue(s).

The legend on the right indicates that replicates of each genotype share the same color, while shape indicates tissue.

Differential expression among tissues and genotypes

We detected 116–2,589 (mean 961–1,398) genes as differentially expressed between genotypes within each tissue (Fig 2A). However, gene expression divergence (i.e. the cumulative differences in expression across all analyzed genes, as measured by 1 –Spearman’s rho, ρ) between genotypes within each tissue was not significantly different among tissues (t-test; Bonferroni-corrected P > 0.8 for all; S1A Fig). Expression divergence tended to be lower between strains derived from the same population (i.e. Swedish strains were more similar to each other than to the Zambian strains and vice versa), although in the hindgut and Malpighian tubule, SU58 was equally or more similar to one or both Zambian strains than to the SU26 strain (Fig 2A). This pattern was not evident in the Malpighian tubule when all genes that could be analyzed in this tissue were included in the analysis (S2 Fig). When we compared expression within the same genotype among tissues, we detected 4,524–5,139 (mean 4,844) genes as differentially expressed between any two tissues (Fig 2B), 50–58% of which were differentially expressed in all pairwise tissue comparisons within a genotype and 1,619–1,880 of which were shared among at least two genotypes, with 1,045 genes differentially expressed among all tissues within all genotypes (S1 Table). Of these shared differentially expressed genes, 1,243–1,594 were consistently upregulated in the same tissue within the same genotype, with 700–924 genes consistently upregulated in the same tissue in all genotypes (S1 Table). Interestingly, overall gene expression divergence within the same genotype between the midgut and Malpighian tubule was significantly lower than gene expression divergence between either of these two tissues and the hindgut (t-test; Bonferroni-corrected P < 5 x 10−5 for both; S1C Fig), suggesting that among these three tissues, expression within the same genetic background is most similar between the Malpighian tubule and the midgut. When we compared gene expression divergence among genotypes within tissues to gene expression divergence within the same genotype among tissues, gene expression divergence was higher among than within tissues (Bonferroni-corrected P = 8.58 x 10−15; Figs 2C and S1A). Thus, expression diverges more within a genotype among tissues than among genotypes within a tissue, suggesting that tissue is more predictive of gene expression than genotype.

10.1371/journal.pgen.1011257.g002 Fig 2 Gene expression divergence among genotypes and tissues.

A) The numbers of differentially expressed (DE) genes between genotypes within the midgut (triangles), hindgut (circles), and Malpighian tubule (squares) are shown above the diagonal, while expression divergence (as measured by 1 –ρ) between genotypes is shown below the diagonal. B) The numbers of differentially expressed genes between the same genotype among midgut (MG), hindgut (HG), and Malpighian tubule (MT) tissues are shown. Dashes indicate missing data. C) Expression divergence among genotypes within the same tissue (across genotype) versus expression divergence between the same genotype among tissues (across tissue). Significance was assessed with a t-test. *** Bonferroni-corrected P < 10−14.

Mode of expression inheritance is highly tissue- and genetic background-specific

In order to understand how the mode of expression inheritance varies among genotypes and tissues, we categorized genes according to their expression in the two parental strains and the respective F1 hybrid into the following categories (see Methods for more details): similar, P1 dominant, P2 dominant, additive, overdominant, and underdominant, with the Swedish strains being P1 and the Zambian strains P2 (Fig 3; S2 and S3 Tables). For all backgrounds and tissues, the similar category (i.e. genes with similar expression in parents and hybrids) was the largest (Fig 3A) and showed the greatest overlap among tissues (S3 and S4 Figs). The basic expression inheritance categories (those with genes showing additive or P1 or P2 dominant expression in hybrids in comparison to parents) were the next largest categories (Figs 3A, S3, and S4), and typically similar numbers of genes were classified into these categories. However, there were some exceptions depending on category, background, and tissue (Fig 3A). For example, 1.5-fold more genes were categorized as dominant in either Swedish strain in comparison to the ZI418 strain in the midgut in comparison to the other tissues, while 3.3–4.6-fold more genes were categorized as dominant in the ZI418 strain in comparison to the SU26 strain in the midgut and Malpighian tubule in comparison to the hindgut. Similarly, 1.8–2.5-fold more genes were categorized as dominant in the ZI197 strain in comparison to either Swedish strain in the midgut than in the hindgut. Genes in these basic inheritance categories were often unique to both the tissue and category (Figs 3B, 3C, S3 and S4), with little overlap within each category across all three examined tissues (Figs 3 and S3). Unsurprisingly, in background combinations for which we only had data for two tissues, the overlap we detected within categories across tissues was higher (S4 Fig). The smallest number of genes were categorized into misexpression categories, i.e. genes showing either over- or underdominance in the hybrid in comparison to the parents (Figs 3A, 3C, S3 and S4). Genes in misexpression categories tended to be tissue-specific with little or no overlap among the examined tissues (Figs 3A, 3C, S3 and S4). Similar to what we observed for basic inheritance categories, certain combinations of genetic backgrounds and tissues showed larger numbers of misexpressed genes than others (Figs 3A, S3, and S4). For example, we detected relatively high levels of misexpression in the SU26xZI418 background in the midgut and Malpighian tubule, but not the hindgut (Fig 3A). Taken together, our results suggest that the mode of expression inheritance is both tissue- and genetic background-specific.

10.1371/journal.pgen.1011257.g003 Fig 3 Dominance divergence and the mode of expression inheritance in examined tissues.

A) The number of genes in each mode of expression inheritance category within the midgut (triangles), hindgut (circles), and Malpighian tubule (squares) at a 1.25-fold change cut-off. Results using alternative cut-offs or for individual tissue analyses can be found in S2 and S3 Tables (see Methods). Dashes indicate missing data. B, C) Upset plots showing unique and overlapping genes within each tissue in the SU26xZI418 background as an example. Upset plots for the other genotypes can be found in S3 and S4 Figs. Horizontal bars represent the total number (num.) of genes in a tissue and inheritance category combination. Vertical bars represent the number of genes in an intersection (intersect.) class. A filled circle underneath a vertical bar indicates that a tissue and inheritance category combination is included in an intersection class. A single filled circle represents an intersection class containing only genes unique to a single tissue and inheritance category combination, while filled circles connected by a line indicate that multiple tissue and inheritance category combinations are included in an intersection class. Genes categorized into B) basic expression inheritance (inherit.), i.e. P1 dominant (P1 dom.), P2 dominant (P2 dom.), and additive (add.) and C) misexpression (misexpress.) categories are shown. Only intersection classes comprised of either a single tissue and inheritance category combination or an inheritance category in all examined tissues are shown. Additional intersection classes and upset plots for genes categorized into the similar category are shown in S3 Fig D) Phenotypic dominance (h) divergence (as measured by 1 –ρ) among backgrounds within the same tissue (across genotype) versus dominance divergence between the same background among tissues (across tissue). Significance was assessed with a t-test. The Bonferroni-corrected P value is shown.

Phenotypic dominance and the mode of expression inheritance

In order to better understand potential variation in the magnitude of phenotypic dominance during expression inheritance, we calculated the degree of dominance, h. In order to compare the magnitude of dominance regardless of which allele was dominant, we calculated h such that values between 0 and 1 or 0 and -1 represent varying degrees of additivity and dominance, with values closer to -1 representing complete dominance of the Swedish background and 1 representing complete dominance of the Zambian background, while values outside this range represent cases of overdominance of the respective background (see Methods for more details). For all genetic backgrounds and tissues, we did not detect any significant difference in the overall magnitude of phenotypic dominance between the two parental backgrounds (t-test, Bonferroni-corrected P > 0.6 for all). For the majority of tissues and genetic backgrounds, we did not detect differences in the magnitude of dominance within the same genetic background between tissues (Bonferroni-corrected P > 0.26 for all comparisons). We only detected a significant difference in the overall magnitude of dominance within the SU26xZI418 background between the midgut and the Malpighian tubule (Bonferroni-corrected P = 0.015), which may be driven by the large amount of misexpression that we detected in this background, particularly in the Malpighian tubule (Fig 3A). Overall, these results suggest that the differences we detected in the mode of expression inheritance among genetic backgrounds and tissues occur on the individual gene level rather than being driven by general, genome-wide changes in dominance. Overall dominance divergence among genetic backgrounds (i.e. the cumulative differences in dominance across all analyzed genes, as measured by 1 –ρ) was not significantly different between the midgut and hindgut (Bonferroni-corrected P = 0.264, S1B Fig), but could not be compared to the Malpighian tubule for which only 2 background combinations were available. When we compared overall dominance divergence among genetic backgrounds within tissues to dominance divergence within the same genetic background among tissues, dominance divergence was significantly higher among than within tissues (Bonferroni-corrected P = 0.012, Fig 3D). We observed a similar pattern when we examined gene expression divergence (Fig 2C), suggesting that in general divergence is higher among tissues than among different genetic backgrounds within a tissue. However, divergence was higher for dominance than for gene expression (t-test, Bonferroni-corrected P < 10−14), suggesting that phenotypic dominance of expression is much less conserved among tissues and genotypes than expression itself, although it is possible that this difference can be explained in part by differences in how each trait was measured as expression was measured in a single genotype but dominance was calculated based on three genotypes.

Genetic basis of expression variation is highly tissue- and genetic background-specific

In order to identify genes in our dataset with any level of cis-regulatory divergence between the parental alleles in any genetic background and tissue, we tested for ASE in genes for which we could distinguish between the parental alleles in the hybrid (see Methods). Of the 4,305–4,592 genes we were able to analyze in all tissues of a genetic background, we detected 80–370 genes showing significant ASE (FDR <5%) depending on genetic background and tissue (Table 1), with a total of 356, 408, 460, and 256 non-redundant genes detected as having ASE in any tissue in the SU26xZI418, SU58xZI418, SU26xZI197, and SU58xZI197 backgrounds, respectively, and a total of 958 genes in all tissues and backgrounds. Within each genetic background 55–86% of genes showing ASE in a particular tissue were unique to that tissue, while, within each tissue, 55–76% of ASE genes were unique to a single genetic background. Indeed, within each genetic background, only 9–73 genes were detected as having ASE in all examined tissues, with backgrounds in which only 2 tissues were examined sharing more ASE genes (Table 1). Thus, allele-specific expression is largely tissue- and genetic background-specific.

10.1371/journal.pgen.1011257.t001 Table 1 ASE genes.

	SU58 vs ZI418 a	SU26 vs ZI418 a	SU58 vs ZI197 a	SU26 vs ZI197 a	
Tissueb	DEP	DEH	ASE	DEP	DEH	ASE	DEP	DEH	ASE	DEP	DEH	ASE	
HG	814	391	252	766	390	235	896	200	172	891	576	370	
MG	491	243	145	406	219	112	685	180	134	545	261	163	
MT	570	114	99	501	108	80	–	–	–	–	–	–	
All	77	20	9	51	24	11	334	69	50	270	133	73	
a A total of 4,035, 4,172, 4,305, and 4,592 genes could be analyzed in the SU58xZI418, SU26xZI418, SU58xZI197, and SU26xZI197 backgrounds, respectively. Results for individual tissue analyses can be found in S4 Table.

b Number of differentially expressed (DE) genes between the parental strains (P) and alleles within the F1 hybrid (H) as well as allele-specific genes (ASE) are shown for hindgut (HG), midgut (MG), Malpighian tubule (MT), and shared across all tissues (All). Dashes indicate missing data.

In order to further understand how the genetic basis of expression varies among genetic backgrounds and tissues, we classified genes in each genetic background and tissue combination into six regulatory categories: “conserved”, “cis-only”, “trans-only”, “cis + trans“, “cis x trans“, “compensatory“, and “ambiguous”[5] (see Methods for more details). The proportion of genes falling into each regulatory category was dependent upon tissue and genetic background, although, in general, when considering genes with non-ambiguous regulatory divergence in all tissues and genetic backgrounds, the largest proportion fell into the trans-only category which contained 2.9–30.6-fold more genes than the cis-only category (Fig 4A). The midgut had a higher proportion of ambiguous genes and a smaller proportion of conserved genes than the other examined tissues (Fig 4A). We detected the most cis-only genes in the hindgut, with 2.3–4.2 fold more genes categorized as cis-only in comparison to the other examined tissues (Fig 4A). In comparison to other genetic backgrounds, the SU58xZI197 background had a higher proportion of ambiguous genes and a smaller proportion of conserved genes as well as 2.2–4.7- and 2.4–10.4-fold fewer genes categorized as cis-only and compensatory, respectively (Fig 4A). Within each genotype, genes with non-ambiguous regulatory divergence were often unique to both the tissue and regulatory category (Figs 4, S5, and S6), with little overlap within each category across all three examined tissues (Figs 4 and S5). Similarly, within each tissue, genes with non-ambiguous regulatory divergence were often unique to a genetic background, with 35–91% of genes unique to a single genetic background within a regulatory class and tissue, while 31–87% of genes were unique to the genetic background and tissue within a regulatory class (S7 Fig). Overlap among all genotypes within a tissue was highest for genes in the trans-only category with 2.8–30.7-fold more shared genes categorized as trans-only in comparison to other non-ambiguous regulatory divergence categories (S7 Fig). Overall, our results suggest that the genetic basis of expression inheritance is both tissue- and genetic background-specific.

10.1371/journal.pgen.1011257.g004 Fig 4 Dominance and the genetic basis of expression variation.

A) The number of genes in each regulatory class within the midgut (triangles), hindgut (circles), and Malpighian tubule (squares). Results for individual tissue analyses can be found in S5 Table (see Methods). Dashes indicate missing data. B) Upset plot showing unique and overlapping genes with non-ambiguous regulatory divergence within each tissue in the SU26xZI418 background. Upset plots for other genotypes can be found in S5 and S6 Figs. Only intersection classes comprised of either a single tissue and regulatory category combination or a regulatory category in all examined tissues are shown. Additional intersection classes are shown in S5 Fig. C) Dominance and D) magnitude of dominance h for genes categorized as cis-only (c, light) and trans-only (t, dark) in each background and tissue. Magnitude of dominance was calculated as the absolute value of dominance h. Only h values with magnitudes of 5 or below are shown. Boxplots including more extreme h values can be found in S8 Fig. Significance was assessed with a t-test. Bonferroni-corrected P values are shown in grey. *** P < 0.005, ** P < 0.01, * P < 0.05.

Phenotypic dominance and the genetic basis of expression variation

Previous studies have found that cis-regulatory variation tends to be more additive [20,27–28], while trans variation tends to be more dominant [27]. In order to better understand the relationship between the genetic underpinnings of expression variation and dominance, we examined phenotypic dominance (h; see Methods) in genes categorized as cis-only or trans-only. Overall dominance in genes categorized as trans-only was often biased towards one parent, with 5 out of 10 tissue and genetic background combinations significantly more dominant in one parental background than the other (Figs 4C and S8; t-test), and the more dominant parent dependent on the tissue and genetic background (Fig 4C). On the other hand, overall phenotypic dominance in genes categorized as cis-only was not significantly biased towards one parental background for any of the examined tissue and genetic background combinations (Fig 4C; t-test, P > 0.5 for all). When we compared the overall magnitude of dominance (as measured by the absolute value of dominance h), trans-only genes were only sometimes more dominant than cis-only genes and this was significant after multiple test correction only in the midgut of the SU26xZI418 background (Fig 4D; t-test, P < 0.05). However, this lack of significance might be due to lack of power, particularly for cis-only genes, of which we detected fewer (Fig 3A). When we included all genes that could be analyzed in each individual tissue and genotype in the analysis, we were able to examine phenotypic dominance in 1.9–3.7-fold more cis-only and 1.6–2.5-fold more trans-only genes (S6 Table). The results, however, remained similar, with no increased detection of differences in dominance (S6 Table), which suggests that although we cannot completely rule it out, these results are unlikely to be due to a lack of statistical power. Thus, although we detected trans-regulatory variants as more dominant and cis-regulatory variants as more additive, which has been reported by previous studies [20, 27–28], we only detected this trend in a single background and tissue. Moreover, the phenotypic dominance of trans- but not cis-regulatory variants tended to be biased toward one parental background.

Functional classification of genes displaying ASE

In order to understand the types of genes showing ASE in our dataset, we tested for an enrichment of gene ontology (GO) biological process and molecular function terms for genes with ASE in each background and tissue. The most commonly enriched GO terms across all backgrounds and tissues were related to oxidoreductase activity (S7 Table). Indeed, we detected at least one oxidoreductase activity term for every genetic background and tissue combination in which we detected enriched GO terms. For genes displaying ASE in all examined tissues within a genotype, we could only detect two enriched GO terms, oxidoreductase activity and response to toxic substance, in the SU26xZI197 background (S7 Table). Thus, despite ASE genes tending to be tissue- and genetic background-specific, in general ASE genes tended to be enriched for genes predicted to be involved in oxidoreductase activity.

ASE genes are enriched for sex-biased gene expression in a context-dependent manner

A previous study on ASE in D. melanogaster found differences in the relative contribution of cis-regulatory effects among genes with different levels of sex bias and among two tissues/body parts [28], while another using hybrids between D. simulans and D. mauritiana found that sexually dimorphic regulatory effects are often in cis [27]. In order to better understand the relationship between sex-biased expression and cis-regulatory variation in the midgut, hindgut, and Malpighian tubule, we categorized genes according to their level of sex bias using data from FlyAtlas2 [41] (see Methods). When considering genes displaying ASE (Table 1), in the Malpighian tubule we detected a significant enrichment of sex-biased genes for both genetic backgrounds (χ2 test, P = 0.001 for both; S8 Table), and in the hindgut we detected an enrichment of sex-biased genes in all genetic backgrounds (P < 0.02 for all; S8 Table) except SU26xZI418 (P = 0.2929; S8 Table). In the midgut, we detected a significant enrichment of sex-biased genes only in the SU26xZI197 background (P = 0.0116; S8 Table), despite the midgut having more sex-biased genes, particularly male-biased genes, than the other tissues (S8 Table). In the tissues and genetic backgrounds where we detected an enrichment of sex-biased ASE genes, the enrichment did not appear to be driven by bias towards one particular sex (<1.63-fold difference in the prevalence of male- and female-biased genes for all), with the exception of the SU26xZI197 background in the midgut, where male-biased genes were 6-fold more prevalent than female-biased genes (S8 Table). It should be noted, however, that only females were used in our experiments to measure ASE. Overall, genes displaying ASE were enriched for sex-biased gene expression, however this enrichment was dependent upon the tissue and genetic background in which they were detected.

The effects of inversions on regulatory variation

Chromosomal inversions have been shown to affect expression in Drosophila [42]. Two inversions that are at high frequency worldwide and in sub-Saharan Africa, respectively [43], were present in our study: In(2L)t in SU26 [19] and In(3R)K in ZI197 [44]. A previous study utilizing the current Malpighian tubule data found that while In(2L)t made a minor contribution to gene expression variation, its presence could not explain the patterns of gene expression detected in F1 hybrids with SU26 as a parent [19]. To better understand how the presence of these inversions potentially affects expression and regulatory variation, we tested for a significant over- or underrepresentation of genes differentially expressed or displaying ASE between parental strains among genes located within these inversions. Genes differentially expressed between parental strains were significantly enriched for genes located within both inversions for all comparisons, including comparisons between strains that did not contain the focal or any inversion (χ2 test, P < 10−15 for all; S9 Table). This finding suggests that the 2L and 3R chromosome arms are enriched for differentially expressed genes, rather than the inversions themselves, which is in line with a previous study that found that linked genetic variation within inversions is what drives differential expression rather than the structural variation itself [45]. For genes displaying ASE (Table 1), we did not find any dearth or enrichment for genes located within either inversion (P > 0.29; S9 Table), with the exception of the SU58xZI418 genetic background in the Malpighian tubule, which was significantly enriched for genes located on the In(2L)t (P = 0.003; S9 Table), despite having the standard chromosomal arrangement. Thus, the presence of the In(2L)t and In(3R)K inversions does not appear to explain the regulatory patterns that we detected.

Tissue specificity varies depending on regulatory type and genetic background

For genetic background combinations for which we had transcriptome data in all three tissues (SU26xZI418 and SU58xZI418), we were able to examine the relationship between regulatory variation and tissue specificity. To do so, for every gene in each strain we calculated the tissue specificity index τ, which ranges in value from 0 to 1, with higher numbers indicating higher tissue specificity. When we compared overall τ among all genetic backgrounds, tissue specificity in ZI418 was higher than in either Swedish background as well as both F1 hybrids, although this difference was not statistically significant for SU26 (Fig 5A). On the other hand, tissue specificity in the Swedish strains was not significantly different from each other or their respective hybrid (Fig 5A). Thus, tissue specificity in F1 hybrids was more similar to the Swedish than the Zambian parent. In order to better understand how tissue specificity varies based on regulatory variation type, we performed pairwise comparisons of τ between genes with trans-only variation, cis-only variation, and conserved gene regulation. Both cis- and trans-regulated genes were significantly more tissue-specific than conserved genes for all strains in all backgrounds (Fig 5B). Interestingly, genes with trans-only regulatory variation were more tissue-specific than genes with cis-only regulatory variation, although this difference was not significant for SU58 and the F1 hybrid in the SU58xZI418 background (P = 0.099 and 0.076, respectively; Fig 5B, S10 Table). Thus, trans effects were more tissue-specific than cis effects in our dataset. Unlike for the genetic basis of expression variation itself (Fig 4), we detected very few tissue-specific or tissue-by-regulatory type interaction effects on tissue specificity (S10 Table). Thus, the genetic basis of regulation (i.e. cis versus trans) and, to a lesser degree, the genetic background are predictive of tissue specificity, while the tissue in which the regulatory variation was detected tends not to be. Indeed, the type of regulatory variation appears to have the largest influence on tissue specificity, as we were able to detect consistent patterns across tissues and genetic backgrounds (Fig 5B).

10.1371/journal.pgen.1011257.g005 Fig 5 Tissue specificity in the SU26xZI418 and SU58xZI418 backgrounds.

A) Overall tissue specificity as measured by τ in the examined strains. Significance was assessed with a t-test and Bonferroni-corrected P values are shown. B) Tissue specificity τ in each strain for genes categorized into cis-only (light), trans-only (tran, dark), and conserved (cons, grey) regulatory (reg.) classes for each genetic background in the midgut (triangles), hindgut (circles), and Malpighian tubule (squares). Significance was assessed with an ANOVA (S10 Table). *** P < 0.005, ** P < 0.01, * P < 0.05, ns not significant (P > 0.05).

Microbiome composition varies depending on tissue and genetic background

Bacterial community composition has been shown to affect host gene expression in the digestive tract depending upon host genotype [46]. In order to better understand the relationship between genetic background and microbiome composition, we performed microbiome sequencing in the midgut and hindgut for the same RNA samples for we which we performed mRNA-seq (see Methods). For all samples, Wolbachia was highly predominant in the bacterial community (10.36–99.37%; S9 Fig). In order to ensure its presence did not mask more subtle differences in diversity, we focus on analyses with Wolbachia removed (Figs 6 and 7), but results including Wolbachia were qualitatively similar (S9 and S10 Figs, S11–S13 Tables), and we did not detect any pattern of relative Wolbachia abundance among genotypes (lmer, P = 0.246; S14 Table). After removal of Wolbachia, Acetobacter, one of the most common D. melanogaster gut microbial taxa [47–48], remained predominant in the midgut (54.7–99.8%; Fig 6). In the hindgut, where the microbiome composition was more diverse, Acetobacter was only dominant in a subset of individuals (1.44–78.52%; Fig 6). Interestingly, we did not detect Lactobacillus, another of the most common gut microbial taxa [47–48]. However, because we performed amplicon sequencing using RNA rather than DNA as the starting material (see Methods), the microbiome composition we detected is representative of metabolic activity rather than presence. Thus, it may be that Lactobacillus was present but its metabolic activity was not high enough for us to detect.

10.1371/journal.pgen.1011257.g006 Fig 6 Composition of the bacterial communities in the midgut and hindgut of each genotype.

Colored sections of each bar show bacterial genera (excluding Wolbachia) with a relative abundance above 5% in each sample. The remaining genera are compiled in the “Others” category. Bacterial community composition including Wolbachia can be found in S9 Fig.

In order to detect differences in gut bacterial community composition, we computed the Bray-Curtis index (see Methods). We detected significant differences in gut bacterial communities between the midgut and hindgut (Figs 6 and 7, S11 Table; PERMANOVA, P = 0.001), which is unsurprising given that these gut regions differ in their pH and associated digestive functions [37]. We also detected a significant effect of the genetic background (i.e. strain) as well as a significant interaction effect between the examined strain and gut region on the bacterial community (Fig 7, S11 Table; PERMANOVA, P ≤ 0.015 for both), suggesting that genetic background affects microbiome composition, and this effect is at least partially tissue-dependent. Indeed, when we examined community composition within each tissue individually, the genetic background significantly influenced the gut bacterial community in the hindgut while it had no significant effect on the structure of the community in the midgut (Fig 7, S11 Table; PERMANOVA, P = 0.002 and P = 0.89 respectively). We also detected significant tissue and genetic background effects on bacterial alpha-diversity (S11 Fig, S12 and S13 Tables; LMM, P = 0.017 and P < 0.001), with the hindgut being more diverse and displaying stronger differentiation between parental and F1 strains (S11 Fig). Thus, the diverse bacterial community of the hindgut offered more possibility for differentiation while the midgut community was dominated by Acetobacter among all samples. In contrast, gene expression in the hindgut was less differentiated among genetic backgrounds but more differentiated from the other tissues while the midgut showed the opposite pattern (Figs 1 and S1C), suggesting that the expression and regulatory variation we detected in these tissues is unlikely to be driven by bacterial community composition. Thus overall, tissue type (i.e. gut region) had the largest impact on microbial community composition and diversity, with genetic background also affecting microbiome variation to a lesser degree, especially within the hindgut; however, these genetic background effects do not appear to be related to host gene expression variation.

10.1371/journal.pgen.1011257.g007 Fig 7 Principal coordinate analysis of bacterial communities in A) both midgut and hindgut samples, B) hindgut, and C) midgut. The legend indicates that replicates of each genotype share the same color, while shape indicates tissue. Wolbachia was excluded from the analysis. Results including Wolbachia can be found in S10 Fig.

Discussion

Using transcriptome data from parental and F1 hybrid D. melanogaster strains from an ancestral and a derived population in the midgut, hindgut and Malpighian tubule, we found that both the genetic basis of expression variation (i.e. cis versus trans) and the mode of expression inheritance (i.e. dominant versus additive) were highly tissue- and genetic background-specific (Figs 3 and 4, Table 1). Previous studies using F1 hybrids in Drosophila have found that genetic background [19,21,28] and body part or tissue [21–22,28] can have large effects on regulatory architecture; however, to our knowledge, this is the first study examining highly spatially and functionally proximate tissues that not only communicate with each other but also functionally and physically interact. Thus, our results demonstrate that even functionally related, interconnected tissues can show highly divergent regulatory architecture among tissues and genetic backgrounds. Indeed, overall gene expression was most similar between the Malpighian tubule and the midgut (Fig 1), despite these tissues being part of the excretory and digestive system, respectively, while the two gut tissues are part of the same alimentary canal. Thus, our results suggest that the level of functional and physical interconnectivity between tissues may not necessarily be predictive of similarity in gene expression or regulatory architecture. Consistent with this interpretation, we detected similar fold-size differences in the proportion of cis-regulated genes in the hindgut versus the Malpighian tubule or midgut (Fig 4A) as has previously been reported in the testes versus the head or ovaries [28]. However, we should note that the tissues we examined in this study, and the midgut in particular, are known to be regionalized [37–38,40] with an estimated 22, 4, and 5 distinct cell types currently described in the midgut, hindgut, and Malpighian tubules, respectively [49–50]; therefore, it is possible that we may have missed some of the more subtle differences in gene regulation that occur among individual regions or cell types.

It has long been thought that regulatory changes and particularly cis-regulatory changes are important during adaptation as they can fine-tune gene expression both temporally and spatially [51–52]. Indeed, we found that genes with trans and cis effects were more tissue-specific than genes with no differential expression regulation (Fig 5), suggesting that regulatory changes between diverged populations are often tissue-specific, which is likely driven by spatial fine-tuning of gene expression. Interestingly and somewhat surprisingly, we found that trans effects were more tissue-specific than cis-effects and this finding was consistent across tissues and genetic backgrounds (Fig 5B). Thus, our results reveal that trans-regulatory changes can be as or, potentially, even more tissue-specific than cis-regulatory changes that occur as populations diverge. In contrast to our findings, a recent study on ASE in two mouse tissues found that tissue-specific genes were largely cis-regulated during population divergence [25]. Indeed, cis-regulatory changes have long been thought to be more common during adaptation due to lower pleiotropy [53]. One of the disadvantages of our methods is that we were unable to assign regulatory effects to their underlying genetic variants and do not know the location or identity of the genetic variants underlying the detected tissue-specific trans effects. Thus, it is possible that the tissue-specific trans effects we detected are driven by cis-regulatory changes in the transcription factors or other regulators driving these trans effects.

Previous studies found that cis-regulatory variation tends to be more additive [20,27,28] than trans variation, which tends to be more dominant [27]. However, we found little evidence for this pattern in our dataset (Fig 4D). The discrepancy between the current study and previous ones may be due to differences in methods or the examined genetic background and/or body parts/tissues, suggesting that differences in additivity and dominance between cis and trans variation may be context-specific. Interestingly, we found that the phenotypic dominance of trans- but not cis-regulatory variation tended to be biased toward one parental background, with the direction of the bias variable among tissues and genetic backgrounds (Fig 4C). The context-dependent nature of this finding suggests that this bias may be driven by one or several trans factors affecting the expression of multiple genes in individual tissues and genetic backgrounds, which underscores the importance of taking genetic background and tissue into account when attempting to identify general patterns and trends in gene expression and its regulation. When we examined divergence in gene expression and phenotypic dominance (i.e. the cumulative differences in each trait across all analyzed genes), we found that divergence was higher among than within tissues (Figs 2C and 3D), suggesting that although both are pervasive, tissue-specific effects outnumber or are larger than genetic background-specific effects, and these effects may be magnified when considering the phenotypic dominance of gene expression, as our findings suggest that it is much less conserved than expression itself.

A previous study in D. melanogaster larvae found that overall developmental (i.e. temporal) gene expression specificity increased during adaptation in a derived population [36]. In contrast, in our dataset the ancestral ZI418 genetic background showed the highest tissue (i.e. spatial) gene expression specificity (Fig 5A); however, it is possible that overall changes in gene expression specificity driven by adaptation are only detectable at the population rather than the individual level. Because we identified regulatory variation between an ancestral and a derived D. melanogaster population that had to adapt to new habitats, some, although not necessarily all, of the regulatory variation we identified may be adaptive. Indeed, one recent study examining ASE between warm- and cold-adapted mouse strains found signs of selection on ASE genes in the cold-adapted mice [25]. Genes we identified as showing ASE included several for which adaptive cis-regulatory divergence has previously been documented, such as MtnA [54], Cyp6g1 [55], Cyp6a20 [19], and Cyp12a4 [19]. We also detected an enrichment of oxidoreductase activity and response to toxic substance among ASE genes (S7 Table), suggesting any genes with adaptive cis-regulatory variation that we detected may be related to these processes. Indeed, the detected selection on Cyp6g1 expression is thought to have been driven by resistance to the pesticide DDT [55], while selection on MtnA is thought to be driven by increased oxidative stress resistance [34,54]. Indeed, the digestive system’s direct interaction with the external environment [37] and the excretory system’s role in detoxification and waste excretion [38] suggest that many of the ASE genes we identified may be candidates for adaptation.

Similar to our findings for the genetic basis of expression variation, the mode of expression inheritance, and phenotypic dominance (Figs 3–5), we detected significant effects of tissue and genetic background on bacterial community composition in our microbiome analysis, although the detected genetic background effects did not appear to explain host gene expression variation (Figs 6, 7, and S11). The endosymbiont Wolbachia, which is known to affect microbiome composition but is not present in the gut lumen [56], was predominant in all of our samples (S9 Fig) but we did not detect the very common Lactobacillus (Fig 6), which suggests that physical abundance within the gut may not be predictive of metabolic activity levels and some bacterial community members may be more or less active than predicted by their physical abundance. Bacterial community composition was highly divergent between the two gut tissues, and the effect of genetic background appeared to be driven by higher diversity in the hindgut, which also showed more differentiation among strains (Figs 6, 7 and S11). Because all flies were reared in the same lab environment, a large portion of the detected bacterial community was likely acquired in the lab. Rearing environment greatly influences bacterial community composition, with communities of lab-reared strains less diverse than their natural-reared counterparts [57]. Thus, it is difficult to draw conclusions about how the genetic background effects we detected might influence bacterial community composition in nature, although genetic differentiation among natural D. melanogaster populations is known to shape bacterial community structure [58].

Overall, our findings yield insight into the evolution of regulatory architecture, the effects of regulatory variation on tissue specificity, the effects of genetic background on expression and microbiome variation, as well as the importance of accounting for context-specificity in evolutionary studies.

Materials and methods

D. melanogaster samples and sequencing

All D. melanogaster strains were reared on cornmeal-molasses medium under standard lab conditions (21°C, 14 hour light: 10 hour dark cycle). mRNA-seq and microbiome sequencing were performed for four isofemale strains, two from Umeå, Sweden (SU26 and SU58) [29] and two from Siavonga, Zambia (ZI418 and ZI197) [30] as well as F1 hybrids between the Swedish and Zambian parental lines (SU58xZI418, SU58xZI197, SU26xZI418, SU26xZI197). The SU58 and ZI418 strains have the standard arrangement for all known chromosomal inversion polymorphisms [19, 44], while SU26 and ZI197 have the standard arrangement with the exception of In(2L)t, which was present in SU26 [19] and In(3R)K, which was present in ZI197 [44]. To determine if inversion status affected our findings, we tested for a significant over- or underrepresentation of genes differentially expressed or displaying ASE between parental strains among genes located within the In(2L)t or In(3R)K inversion using a χ2 test. Reciprocal F1 hybrids were generated in both directions (i.e. parental genotypes were switched) by crossing 2–3 virgin females of one line with 4–5 males of the other line. Crosses were carried out in 8–13 replicate vials and parental strains were similarly reared (2–3 females and 3–5 males per vial with 8–12 replicate vials) in order to control for rearing density among genotypes.

Midguts (from below the cardia to the midgut/hindgut junction, 20 per biological replicate) and hindguts (from the midgut/hindgut junction to the anus, 60 per biological replicate) were dissected from 6-day-old females in cold 1X PBS and stored in RNA/DNA shield (Zymo Research Europe; Freiburg, Germany) at -80°C until RNA extraction. For F1 hybrids, half of the tissues were dissected from each of two reciprocal crosses in order to avoid potential parent-of-origin effects, although such effects are expected to be absent or very rare in D. melanogaster [6,59]. RNA was extracted from three biological replicates per genotype and tissue type (48 samples in total) with the RNeasy Mini kit (Qiagen; Hilden, Germany) as directed by the manufacturer. mRNA-seq and microbiome sequencing were performed using the same RNA extractions. Poly-A selection, fragmentation, reverse transcription, library construction, and high- throughput sequencing was performed by Novogene (Hong Kong) using the Illumina HiSeq 2500 platform (Illumina; San Diego, CA) with 150-bp paired reads. Malpighian tubule 125-bp paired read data for SU58, SU26, ZI418 and F1 hybrids (SU58xZI418, SU26xZI418), which was composed of 2 biological replicates per genotype (10 in total; 58 libraries in total across all tissues) was downloaded from Gene Expression Omnibus (accession number GSE103645).

Microbiome sequencing and analysis

Reverse transcription was carried out to generate complementary DNA (cDNA) which was used for amplicon sequencing targeting the V4 region of the 16S rRNA bacterial gene. First, template RNA was cleaned of potential residual genomic DNA with the PerfeCta DNase I (Quantabio; Beverly, MA) following the manufacturer’s instructions. Reverse transcription was performed using the FIREScript RT cDNA Synthesis (Solis BioDyne; Tartu, Estonia) with specific bacterial primers, 515F (5’-GTGYCAGCMGCCGCGGTAA-3’) and 806R (5’- GGACTACNVGGGTWTCTAAT-3’), also following the manufacturer’s instructions. The V4 region of the 16S rRNA gene was sequenced from the resulting cDNA on an Illumina Miseq platform using the 515F and 806R primer pair. Using the R package DADA2 (version 1.26.0) [60], Amplicon Sequence Variants (ASVs) were inferred after trimming (length of 240nt for forward reads and 180nt for reverse reads). Dereplication and chimera removal were performed using default parameters of DADA2. Each ASV was assigned taxonomically using the Silva classifier (version 138.1) [61]. ASVs assigned to the Eukaryotic and Archeal kingdoms were removed. Given that the gut bacterial community was highly dominated by one ASV assigned to the genus Wolbachia (S9 Fig), a known intracellular symbiont of Drosophila melanogaster, we chose to remove it for further statistical analysis, revealing the underlying diversity in the gut bacterial community.

All statistical analyses were performed in R-4.2.2 and each graph was generated with the ggplot2 package [62]. The composition of the bacterial gut community was analyzed using the Phyloseq package [63]. ASVs not present in more than 6.25% (3 replicates/48 samples = 0.00625) of the samples were removed for visualization purposes but kept in the data for the remaining analyses. Differences in beta-diversity were tested with permutational multivariate analysis of variance (vegan package version 2.6–4) [64] on a Bray-Curtis dissimilarities matrix and Principal coordinate analyses (PCoA) was performed for visualization using the vegan package. Differences in bacterial alpha diversity (species richness, Shannon index, Simpson index and inverse Simpson index generated with vegan package) were tested with linear mixed models (lmer, lme4 package) [65] and pairwise comparisons were tested following the Tukey method (emmeans package) [66]. The RNA extraction batch had no significant effect on differences in alpha and beta-diversity of the bacterial community.

mRNA-seq analyses

Reference genomes for each parental strain were constructed using published genome sequence assemblies of SU26 and SU58 [29], and ZI197 and ZI418 [44, 67] as described in [19]. Briefly, if a nucleotide sequence difference on the major chromosome arms (X, 2R, 2L, 3R, 3L) occurred between a parental strain and the D. melanogaster reference genome (release 6) [68], the parental nucleotide variant was included in the new reference transcriptome. If the parental sequence contained an uncalled base (“N”), the reference sequence was used. All transcribed regions (including rRNAs, non-coding RNAs, and mRNAs) were then extracted from each parental reference genome using FlyBase annotation version 6.29 [68]. For each parental strain library, mRNA-seq reads were mapped to the corresponding parental reference genome. In order to prevent mapping bias for genes with greater sequence similarity to one of the parental reference genomes, reads for F1 hybrids were mapped to the combined parental reference genomes.

Reads were mapped to the reference transcriptomes using NextGenMap [69] in paired-end mode. Read pairs matching more than one transcript of a gene were randomly assigned to one of the transcripts of that gene. For downstream analyses, we analyzed the sum of read counts across all of a gene’s transcripts (across all annotated exons), i.e. on the individual gene-level. To identify genes with poor mapping quality, for each parental transcriptome, we simulated mRNA-seq data with 200 reads per transcript and either 125 bp or 150 bp reads, then mapped the reads back to the corresponding transcriptome. Genes for which more than 5% of reads mapped incorrectly were removed from the analyses of the corresponding read length (125 bp for Malpighian tubule and 150 bp for midgut and hindgut). Library size ranged from 34.7 to 55.0 million paired end reads, 97.0–98.6% of which could be mapped (S15 Table).

ASE and mode of expression inheritance analyses were performed within individual tissues as well as for all tissues together. Analyses for individual tissues were qualitatively similar to our analyses including all tissues; therefore, we focus in the main text on analyses including all tissues (S1 and S2 Data) and have included individual tissue analyses as Supplementary material (S3–S5 Data, S3–S6 Tables). To standardize statistical power across all libraries included in the analysis, we held the total number of mapped reads constant by setting the maximum number of mapped reads per sample to that of the library with the fewest mapped reads and randomly subsampling reads (without replacement) until the total number of mapped reads for each sample equaled the maximum. The number of reads we subsampled for each dataset were as follows: 34,009,757 in midgut, 31,611,417 in hindgut, and 30,820,759 in Malpighian tubule as well as for analyses including all tissues. We identified differentially expressed genes using a negative binomial test as implemented in DESeq2 [70]. Gene expression divergence between two strains or tissues was calculated as 1 –Spearman’s ρ between the mean normalized gene counts of the two. Significant differences in divergence were assessed with a t-test. To be considered as expressed in our dataset, we required that a gene have a minimum of 15 reads in each sample, which resulted in 7,684, 8,209, 7,675, and 6,894 genes in the midgut, hindgut, Malpighian tubule, and all tissues, respectively, that could be used in analyses.

Calculation of tissue specificity and phenotypic dominance h

We calculated normalized gene counts for each sample using DESeq2 [70] for genes expressed in all tissues. We then used the normalized gene counts to calculate the tissue-specificity index tau, τ, [71] for each genotype for which we had data from all three examined tissues and were able to examine tissue specificity for 3,338 genes expressed in all genetic backgrounds and tissues. The degree of phenotypic dominance (h) was calculated for each set of parental strains and their respective F1 hybrid (4 genetic background combinations in total) as: h=2XF1−XZI−XSUXZI−XSU, (1)

where XZI, XSU, and XF1 represent the mean normalized gene count across all replicates for the Zambian parental strain, the Swedish parental strain, and the F1 hybrid, respectively [72] in each set of background combinations. This equation for phenotypic dominance yields values between -1 (complete dominance of the Swedish background) and 1 (complete dominance of the Zambian background), which allows for a simple and intuitive comparison of the magnitude of dominance between the two backgrounds but differs slightly from how phenotypic dominance is calculated by other methods (for example, see [73]). Divergence in the degree of phenotypic dominance (h) between two genetic backgrounds or tissues was calculated as 1 –Spearman’s ρ between the two. Significant differences in divergence were assessed with a t-test.

Inference of the mode of expression inheritance

To infer the mode of expression inheritance in F1 hybrids, we compared F1 hybrid expression to parental expression and classified genes into six categories: “similar,” “P1 dominant”, “P2 dominant”, “additive,” “overdominant,” and “underdominant” [5]. To do so, we compared the fold-change difference in expression as calculated by DESeq2 [70] for each gene between genotypes to a fold-change cutoff threshold. Genes where all expression differences were below the cutoff were classified as “similar”, while genes for which the expression difference was greater than the cutoff between the hybrid and only one parent were classified as dominant for that parent. Genes were categorized as additive if the expression differences between the hybrid and both parents was above the cutoff and the hybrid expression was between the expression of the two parental strains, or if the difference in expression between the two parents was above the cutoff and hybrid expression was between the two parental strains. Genes were categorized as overdominant if the expression difference between the hybrid and both parents was above the cutoff and hybrid expression was greater than that of both parents. A gene was categorized as underdominant if the expression difference between the hybrid and both parents was above the cutoff and hybrid expression was lower than that of both parents. We employed three fold-change cutoffs (1.25, 1.5, and 2) as well as a negative binomial test [70] and a 5% FDR cutoff, for which we also included an ambiguous category for genes that did not fit into the other categories. In the main text, we focus on the 1.25-fold cutoff as i) the relative proportion of genes falling into each of the non-similar categories was qualitatively similar for all cutoffs (S2 Table), ii) a fold-change cutoff (rather than a statistical test) should avoid bias in detecting differential expression between alleles/genes with higher expression, as the power of statistical tests increases with increasing read counts, iii) the 1.25-fold cutoff has been employed in several previous studies with the justification that most of the significant expression differences detected between samples tend to be of this magnitude [5,20,74], and iv) previous work using the Malpighian tubule data we use here empirically determined it to be a reasonable cutoff for this analysis [19]. The results for the other cut-offs and individual tissues are provided in S2 and S3 Tables.

ASE analysis

In order to detect expression differences between the two alleles in each F1 hybrid, we compiled lists of diagnostic SNPs that could be used to distinguish between transcripts for each pairwise combination of parental alleles in each examined tissue (S16 Table). To do so, we first compared the two parental reference genome sequences over all transcribed regions annotated in the D. melanogaster reference genome (version 6.29) [68] to compile an initial list of diagnostic SNPs. In order to exclude sites with potential residual heterozygosity or sequencing errors, for each tissue, we required that all SNPs inferred from the genome sequences be confirmed in the parental mRNA-seq data with a coverage of ≥ 20 reads and the expected variant in ≥ 95% of the mapped reads of each parent. Next, we called new SNPs from the parental mRNA-seq data if a site was not polymorphic (or contained an N) in the parental genome sequences, but had ≥ 20 mapped reads in each parent with ≥ 95% having the same base in one parent but a different base the other parent. The total number of high-confidence diagnostic SNPs meeting these criteria was 66,030–89,497, 74,388–105,634, and 59,294–60,668 for midgut, hindgut, and Malpighian tubule, covering 6,937–7,474, 7,399–8,166, and 6,394–6,412 genes, respectively.

To assess ASE in the F1 hybrids, we used the mapping data described in the mRNA-seq analyses section above, but used only reads containing at least one diagnostic SNP (i.e. reads that could be assigned to a parental allele). As described above, counts were summed over all transcripts of a gene and the ASE analysis was carried out on a per gene basis. To standardize statistical power between genetic background combinations and/or tissues while maximizing the number of reads that could be included in our analysis, the maximum number of diagnostic reads per sample (i.e. the maximum number of reads for 2 alleles) was set to that of the F1 hybrid with the fewest diagnostic reads. For all other samples, reads were randomly subsampled (without replacement) until the total number of diagnostic reads equaled the maximum for F1 hybrids or half of the maximum for parents. The maximum number of diagnostic reads was set to 15,446,286, 11,058,789, 12,477,714, and 11,058,789 for midgut, hindgut, Malpighian tubule, and all tissues, respectively. We tested for differences in allelic expression using a negative binomial test as implemented in DESeq2 [70], using only genes with a minimum of 15 diagnostic reads for each allele replicate, resulting in a total of 5,060–5,590, 5,650–6,141, 5,097–5,133, and 4,035–4,592 genes depending on genetic background combination that could be analyzed in the midgut, hindgut, Malpighian tubule, and all tissues, respectively, of which 4,228, 4,800, 4,397, and 2,845 genes could be analyzed in all genetic background combinations. In the main text, we focus on the 2,845 genes that could be directly compared across all genetic background combinations and tissues, although results for individual tissues and genetic background combinations were qualitatively similar (S5 Table).

Inference of the genetic basis of expression variation

We determined the genetic basis of expression variation for each gene using the outcome of three statistical tests: a negative binomial test for differential expression between the two parental strains, a negative binomial test for ASE in the F1 hybrid, and a Cochran–Mantel–Haenszel (CMH) test of the ratio of expression between the two parents and the ratio between the two alleles in the hybrid. For all tests, P-values were adjusted for multiple testing [75] and an FDR cutoff of 5% was used to define significant differences. We employed the same subsampling procedure as described in the ASE analysis section above in order to balance statistical power between parents and hybrids. We classified genes into regulatory classes [5] as follows: “conserved” genes showed no significant difference in any test; “all cis” genes showed significant ASE in hybrids and significant DE between parents, but the CMH test was not significant; “all trans” genes showed significant DE between the parents and a significant CMH test, but no ASE; “compensatory” genes had no DE between parents, but showed significant ASE in hybrids and a significant CMH test; “cis + trans” genes were significant result for all three tests with the expression difference between the parents greater than the difference between the two alleles in the hybrid; “cis × trans” genes also had three significant tests, but the expression difference between the parents was less than the difference between the two alleles in the hybrid; and “ambiguous” genes were significant for only one test.

Gene set enrichment analysis

We used InterMine [76] to search for an enrichment of gene ontology (GO) biological process and molecular function terms for genes displaying ASE in each genetic background and tissue as well as in all tissues.

Sex-biased gene analysis

In order to calculate sex bias for each gene in each examined tissue, we downloaded male (M) and female (F) FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values from FlyAtlas2 [41] and calculated sex bias as log2(FPKMM/FPKMF). Based on this log fold-change (LFC), we categorized genes as either sex-biased (LFC ≥ 0.5 or LFC ≤ -0.5) or unbiased (LFC ≤ 0.5 or LFC ≥ -0.5). We tested for a significant over- or underrepresentation of sex-biased genes among genes displaying ASE using a χ2 test with a Benjamini-Hochberg multiple test correction. To better understand how these sex-biased genes were distributed among different levels of sex bias, we further categorized genes as strongly female-biased (FS; LFC ≤ -1.5), moderately female-biased (FB; -0.5 ≤ LFC ≤ -1.5), unbiased (UB; -0.5 > LFC > 0.5), moderately male-biased (MB; 0.5 ≥ LFC ≥ 1.5), or strongly male-biased (MS; LFC ≥ 1.5; S8 Table).

Supporting information

S1 Fig Expression and dominance (h) divergence within and among tissues.

A) Expression and B) dominance (h) divergence among genotypes within the midgut (MG), hindgut (HG), and Malpighian tubule (MT) versus divergence between the same genotype among tissues (across). C) Expression and D) dominance (h) divergence within the same genotype among tissues. A–C) Significance was assessed with a t-test. D) Significance was not assessed due to the low number of comparisons. Bonferroni-corrected P values are shown. * P < 0.05, ** P < 5 x 10−5, *** P < 10−14, ns not significant, nt not tested.

(PDF)

S2 Fig Differential expression and divergence within tissues.

The total number of differentially expressed (DE) genes between genotypes within the A) hindgut (HG), B) midgut (MG), and C) Malpighian tubule (MT) are shown above the diagonal, while expression divergence (as measured by ρ subtracted from one) between genotypes is shown below the diagonal. Analysis was performed in each tissue individually. The numbers of genes that could be included in the analysis for each tissue were 8,209 in the hindgut, 7,684 in the midgut, and 7,675 in the Malpighian tubule.

(PDF)

S3 Fig Mode of expression Inheritance in SU26xZI418 and SU58xZI418 backgrounds.

Upset plots showing unique and overlapping genes within the hindgut (circles), midgut (triangles), and Malpighian tubule (squares) in the A,C,E) SU26xZI418 or B,D,F) SU58xZI418 backgrounds. Horizontal bars represent the total number (num.) of genes in a tissue and inheritance category combination. Vertical bars represent the number of genes in an intersection class. A filled circle underneath a vertical bar indicates that a tissue and inheritance category combination is included in an intersection class. A single filled circle represents an intersection class containing only genes unique to a single tissue and inheritance category combination. Filled circles connected by a line indicate that multiple tissue and inheritance category combinations are included in an intersection class. Genes categorized into A,B) basic expression inheritance (inherit.), i.e. P1 dominant (P1 dom.), P2 dominant (P2 dom.), and additive (add.), C,D) misexpression (misexpress.), and E,F) similar categories are shown.

(PDF)

S4 Fig Mode of expression inheritance in SU26xZI197 and SU58xZI197 backgrounds.

Upset plots showing unique and overlapping genes within the hindgut (circles) and midgut (triangles) in the A,B,E) SU26xZI197 or C,D,F) SU58xZI197 backgrounds. Horizontal bars represent the total number (num.) of genes in a tissue and inheritance category combination. Vertical bars represent the number of genes in an intersection class. A filled circle underneath a vertical bar indicates that a tissue and inheritance category combination is included in an intersection class. A single filled circle represents an intersection class containing only genes unique to a single tissue and inheritance category combination. Filled circles connected by a line indicate that multiple tissue and inheritance category combinations are included in an intersection class. Genes categorized into A,C) similar, B,D) basic expression inheritance (inherit.), i.e. P1 dominant (P1 dom.), P2 dominant (P2 dom.), and additive (add.), and E,F) misexpression (mis-express.) categories are shown.

(PDF)

S5 Fig Genetic basis of expression inheritance in SU26xZI418 and SU58xZI418 backgrounds.

Upset plots showing unique and overlapping genes with non-ambiguous regulatory divergence in the hindgut (circles), midgut (triangles), and Malpighian tubule (squares) in the A) SU26xZI418 or B) SU58xZI418 backgrounds. Horizontal bars represent the total number of genes in a tissue and regulatory category combination. Vertical bars represent the number of genes in an intersection class. A filled circle underneath a vertical bar indicates that a tissue and inheritance category combination is included in an intersection class. A single filled circle represents an intersection class containing only genes unique to a single tissue and regulatory category combination. Filled circles connected by a line indicate that multiple tissue and regulatory category combinations are included in an intersection class.

(PDF)

S6 Fig Genetic basis of expression inheritance in SU26xZI197 and SU58xZI197 backgrounds.

Upset plots showing unique and overlapping genes with non-ambiguous regulatory divergence in the hindgut (circles) and midgut (triangles) in the A) SU26xZI197 or B) SU58xZI197 backgrounds. Horizontal bars represent the total number (num.) of genes in a tissue and regulatory category combination. Vertical bars represent the number of genes in an intersection class. A filled circle underneath a vertical bar indicates that a tissue and inheritance category combination is included in an intersection class. A single filled circle represents an intersection class containing only genes unique to a single tissue and regulatory category combination. Filled circles connected by a line indicate that multiple tissue and regulatory category combinations are included in an intersection class.

(PDF)

S7 Fig Genetic basis of expression inheritance across examined tissues and backgrounds.

Shown are A,B) unique and C,D) overlapping genes in each regulatory category. Shown are A) the number of genes unique to each tissue within each regulatory category and genetic background, B) the number of genes unique to each genetic background and tissue within each regulatory category, C) the number of genes in each regulatory category detected in all examined tissues for each genetic background, and D) the number of genes in each regulatory category detected in all genetic backgrounds for each tissue. Asterisks (*) indicate comparisons using only C) two tissues or D) two genetic backgrounds.

(PDF)

S8 Fig All dominance in cis-only versus trans-only genes.

A) Dominance and B) magnitude of dominance h for genes categorized as cis-only (c, light) and trans-only (t, dark) in each background and tissue. Significance was assessed with a t-test. Bonferroni-corrected P values are shown in grey. *** P < 0.005, ** P < 0.01, * P < 0.05, ms P marginally significant after multiple test correction (P < 0.1), ns P not significant after multiple test correction.

(PDF)

S9 Fig Composition of the bacterial communities in the midgut and hindgut of each genotype, including Wolbachia ASVs.

Colored sections of each bar show bacterial genera with a relative abundance superior to 5% in each sample. The remaining genera are compiled in the “Others” category.

(PDF)

S10 Fig Principal coordinate analysis of bacterial communities in A) both midgut and hindgut samples, B) hindgut, and C) midgut, including Wolbachia ASVs.

The legend indicates that replicates of each genotype share the same color, while shape indicates tissue.

(PDF)

S11 Fig Shannon diversity index of the bacterial community in the midgut and hindgut excluding (A) or including (B) Wolbachia ASVs.

* indicates significant differences of the Shannon index between groups (lmer, P < 0.05).

(PDF)

S1 Table Shared differentially expressed genes among tissues and genotypes.

The numbers of overlapping differentially expressed (DE) genes within or among genotypes and/or tissues for all DE genes (All) or genes upregulated in the midgut (MGup), hindgut (HGup), and Malpighian tubule (MTup) tissues are shown. Only genotypes for which data was available in all examined tissues are shown.

(XLSX)

S2 Table Mode of expression inheritance in combined tissue analysis.

Numbers of genes in each mode of expression inheritance category within the hindgut, midgut, and Malpighian tubule at a 1.25-, 1.5-, and 2-fold change or 5% FDR cut-off (see Methods) are shown. 6,894 genes could be included in the analysis. The ambiguous category is only necessary for the 5% FDR cutoff and comprises genes which could not be assigned into other categories.

(XLSX)

S3 Table Mode of expression inheritance in individual tissue analysis.

Numbers of genes in each mode of expression inheritance category within the hindgut, midgut, and Malpighian tubule at a 1.25-, 1.5-, and 2-fold change or 5% FDR cut-off (see Methods) are shown. 7,684, 8,209, and 7,675 genes could be included in the analysis in the midgut, hindgut and the Malpighian tubule, respectively. The ambiguous category is only necessary for the 5% FDR cutoff and comprises genes which could not be assigned into other categories.

(XLSX)

S4 Table ASE genes identified in individual tissue analyses.

Number of differentially expressed (DE) genes between the parental strains (P) and alleles within the F1 hybrid (H) as well as allele specific genes (ASE) are shown for hindgut (HG), midgut (MG), Malpighian tubule (MT), and shared across all tissues (All). Dashes indicate missing data.

(XLSX)

S5 Table The genetic basis of expression inheritance in individual tissue analyses.

Numbers of genes in each regulatory category within the hindgut, midgut, and Malpighian tubule are shown. 4,228, 4,800, and 4,397 genes could be included in the analysis in the midgut, hindgut and the Malpighian tubule, respectively.

(XLSX)

S6 Table Phenotypic dominance in all cis and all trans groups including all genes that could be analyzed in each individual tissue and genotype.

The mean and mean of the absolute value of phenotypic dominance (h) are shown. Significance was assessed using a t-test. Significant P-values are in bold and values non-significant after multiple test correction are shown in grey.

(XLSX)

S7 Table Enriched GO terms in genes with ASE.

Enriched molecular function and biological process GO terms for genes showing ASE in the Malpighian tubule (MT), midgut (MG), hindgut (HG), and all examined tissues. Number (num) of contributing terms and Holm-Bonferroni-adjusted P-values are shown. Genetic background and tissue combinations with no enriched GO terms are not shown.

(XLSX)

S8 Table Sex bias in genes displaying ASE.

The number of genes categorized as strongly female-biased (FS; LFC ≤ -1.5), moderately female-biased (FB; -0.5 ≤ LFC ≤ -1.5), unbiased (UB; -0.5 > LFC > 0.5), moderately male-biased (MB; 0.5 ≥ LFC ≥ 1.5), strongly male-biased (MS; LFC ≥ 1.5), and generally sex-biased (SB; LFC ≥ 0.5 or LFC ≤ -0.5) are shown. Significant deviations from the expected number of sex-biased genes (SB EXP) were assessed with a χ2 test. P- and Benjamini-Hochberg-corrected P-values are shown.

(XLSX)

S9 Table The effect of In(2L)t and In(3R)K inversion status.

The number of genes differentially expressed (DE) or displaying allele specific expression (ASE) between parental strains along with the number of genes located within an inversion are shown. The strain containing an inversion is shown in parentheses. Significant deviations from the expected number of DE or ASE genes were assessed with a χ2 test.

(XLSX)

S10 Table ANOVA results for pairwise comparisons between genes categorized as cis-only, trans-only, and conserved.

Results for ANOVAs with τ of the indicated strain as the response variable and regulatory (reg.) variant type, tissue, and the interaction between them as factors.

(XLSX)

S11 Table PERMANOVA results on comparison of Bray-Curtis distances of the bacterial community including or excluding Wolbachia ASVs.

(XLSX)

S12 Table Summary statistics of the LMER ANOVA comparing Shannon indices of the bacterial community including or excluding Wolbachia ASVs.

The mixed models included the genotype and tissue as fixed factors and the batch and sample ID as random factors. The interaction term in between the two fixed effect did not have a significant effect on the Shannon indices and was subsequently removed from the models. LMER: Shannon index ~ genotype + tissue + 1| batch + 1| sample ID.

(XLSX)

S13 Table Pairwise comparisons of the Shannon index for the bacterial communities including or excluding Wolbachia ASVs.

The comparisons were performed following the Tukey method and the P-values were adjusted via the Benjamini-Hochberg (BH) method.

(XLSX)

S14 Table Results of lmer comparing the relative abundance of Wolbachia in the bacterial communities.

(XLSX)

S15 Table Library size and mapping efficiency for all mRNA-seq libraries.

To prevent mapping bias, all libraries were simultaneously mapped to at least one pair of parental strains (P1 and P2). The number of mapped, paired mapped, unmapped, and discarded reads as well as total reads, read pairs, and proportion (prop) of mapped reads are shown.

(XLSX)

S16 Table Number of diagnostic SNPs in each tissue and comparison.

Shown are the number of diagnostic SNPs determined from the respective reference genome (Ref) and with the mRNA-seq data included as well as the number of genes covered by the diagnostic SNPs.

(XLSX)

S1 Data ASE gene counts and expression in analyses including all tissues.

(XLSX)

S2 Data Phenotypic dominance, sex-biased gene expression, and τ in analyses including all tissues.

(XLSX)

S3 Data ASE gene counts and expression, and phenotypic dominance in midgut analyses.

(XLSX)

S4 Data ASE gene counts and expression, and phenotypic dominance in hindgut analyses.

(XLSX)

S5 Data Gene counts, expression and phenotypic dominance in Malpighian tubule analyses.

(XLSX)

S6 Data ASV and taxonomy tables of the bacterial communities in the midgut and hindgut.

(XLSX)

We thank Hilde Lainer for excellent technical assistance in the lab as well as Dr. Grit Kunert (Department of Biochemistry, Max-Planck-Institute for Chemical Ecology) for advice on statistical models. We also thank the LMU Evolutionary Biology department for helpful suggestions and discussions.

10.1371/journal.pgen.1011257.r001
Decision Letter 0
Mackay Trudy F.C. Consulting Editor - PLoS Genetics
Williams Scott M. Section Editor
© 2024 Mackay, Williams
2024
Mackay, Williams
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version0
13 Jun 2024

Dear Dr. Glaser-Schmitt, 

Thank you very much for submitting your Research Article entitled 'Pervasive tissue-, genetic background-, and allele-specific gene expression effects in Drosophila melanogaster' to PLOS Genetics.

The manuscript was fully evaluated at the editorial level and by independent peer reviewers. The reviewers appreciated the attention to an important topic but identified some concerns that we ask you address in a revised manuscript.

We therefore ask you to modify the manuscript according to the review recommendations. Your revisions should address the specific points made by each reviewer.

In addition we ask that you:

1) Provide a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript.

2) Upload a Striking Image with a corresponding caption to accompany your manuscript if one is available (either a new image or an existing one from within your manuscript). If this image is judged to be suitable, it may be featured on our website. Images should ideally be high resolution, eye-catching, single panel square images. For examples, please browse our archive. If your image is from someone other than yourself, please ensure that the artist has read and agreed to the terms and conditions of the Creative Commons Attribution License. Note: we cannot publish copyrighted images.

We hope to receive your revised manuscript within the next 30 days. If you anticipate any delay in its return, we would ask you to let us know the expected resubmission date by email to plosgenetics@plos.org.

If present, accompanying reviewer attachments should be included with this email; please notify the journal office if any appear to be missing. They will also be available for download from the link below. You can use this link to log into the system when you are ready to submit a revised version, having first consulted our Submission Checklist.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Please be aware that our data availability policy requires that all numerical data underlying graphs or summary statistics are included with the submission, and you will need to provide this upon resubmission if not already present. In addition, we do not permit the inclusion of phrases such as "data not shown" or "unpublished results" in manuscripts. All points should be backed up by data provided with the submission.

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Please review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript. If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

PLOS has incorporated Similarity Check, powered by iThenticate, into its journal-wide submission system in order to screen submitted content for originality before publication. Each PLOS journal undertakes screening on a proportion of submitted articles. You will be contacted if needed following the screening process.

To resubmit, you will need to go to the link below and 'Revise Submission' in the 'Submissions Needing Revision' folder.

Please let us know if you have any questions while making these revisions.

Yours sincerely,

Trudy F.C. Mackay

Consulting Editor - PLoS Genetics

PLOS Genetics

Scott Williams

Section Editor

PLOS Genetics

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors examine gene expression in three related tissues in individuals from two geographically separated populations plus F1 individuals to study both tissue and allele specific effects on gene expression. The authors find a large number of differentially expressed genes between tissues and backgrounds and that tissue specific effects were larger than background effects. This finding provides insight into the role of cis and trans effects in differential gene expression. The main takeaway is that there is a lot of variability between tissues and background in terms of genes that show ASE and gene regulation. The regulatory variation vs tau comparison is quite interesting.

The methods used in the paper are straightforward and have been used in many other similar analyses.

Comments:

Figure1: the color scheme used in this figure makes it difficult to distinguish some of the genotypes, in particular the two ZI genotypes from the SU58xZI crosses. Please make these colors more distinct.

Figure 2 B: What is the overlap between the gene lists summarized in this table? How much sharing is there between the lists of DE genes between comparisons of the same type with respect to population or the F1s?

The microbiome composition is a nice addition and while I agree that removing Wolbachia seems reasonable in general I have concerns about how this affects the analysis. Is the difference in % Wolbachia significant between either the Sweden and Zambian populations or either parent and the F1 for the midgut?

Number of cell types in each tissue - this is mentioned briefly in the discussion, but I think a bit more context about the tissues would be appropriate here. Additional information about the number of cell types known in each of these tissues would help to give context to the findings in the paper.

The authors point out that there are differences in the inversion status between some of the lines, but do not address if there is any effect of the inversion in the regulatory variation. Is there, for example, an increase in trans effects within the In(2L)t region for comparisons involving SU26?

Reviewer #2: In this study, Glaser-Schmitt et al. examine variation in transcript abundance in three tissues (two functionally related ones and one somewhat less related) and in four genetic backgrounds of D. melanogaster to elucidate the evolution of regulatory architecture. To evaluate potential interactions between gene expression variation and microbiome communities, the authors also characterize bacterial composition in the two gut tissues. In both transcriptome and microbiome data, the authors detected pervasive genetic background- and tissue-specific effects, with the latter being typically stronger than the former. Analyses of allele specific expression (ASE) showed that genes with trans-effects are the most abundant regulatory class among genes with non-ambiguous regulatory variation. Not surprising, tissue specificity was higher for both cis- and trans-regulated genes than for genes with conserved expression across all genetic backgrounds. However, genes with trans-effects were more tissue-specific than genes with cis-effects. This observation suggests that variation in gene expression due to trans-regulatory elements can be spatially fine-tuned as well as (or even better than) that regulated by cis elements.

This is a well-written manuscript based on mRNAseq data analysis. The analyses are sound, and the figures of the manuscript are elaborated and successfully encapsulate the main points of the study in a clear and informative way, which is remarkable when we consider the large amount of new data reported. The results have a large descriptive component and yet they are interesting and exciting. This study fills a current gap in the knowledge of the field and will spark new research on the genetic basis of regulatory evolution. I don’t have any major issues, so I am recommending acceptance for publication in PLoS Genetics with minor revisions. I have just a few minor points and some editorial comments (see below).

1. The term dominance (h) divergence is used to describe variation in degree of dominance among genotypes (within the same tissue) but also among tissues (within the same genotype). It is not clear, however, how dominance divergence is calculated. Furthermore, the word divergence is usually reserved for differences between species but in this study, it also refers to differences among tissues so, in my view, it is a bit confusing. Perhaps the authors would consider alternative terms. “Dominance distance”? Along these lines, what is the difference between “Dominance” and “Magnitude Dominance” (Fig. 4)? Additional clarification is needed.

2. Following the previous point, the terms “gene expression divergence” are used to refer to differences in gene expression among tissues. However, “gene expression divergence” has been extensively used to describe differences in gene expression between species so it would be better to find an alternative. Maybe “gene expression distance?” Not an easy task.

3. This study focuses on gene expression in females. Is there sexually dimorphic gene expression in the 3 tissues examined in D. melanogaster? How much? I wonder about the relative contribution of genes with different degrees of sex-biased expression to the main findings of the study. One could potentially analyze male-, female- and nonsex-biased genes separately, and see to what extent conclusions hold. This may be relevant, considering that the different categories of sex-biased expression genes have been shown to have, for instance, different rates of protein evolution. The downside to this approach is perhaps the potential loss of statistical power. Just a thought for the authors to consider.

4. A technical comment. One of the observations of this study is that divergence for dominance is significantly higher than for gene expression (Bonferroni-corrected P < 10-14), which suggests that phenotypic dominance of expression is much less conserved among tissues and genotypes than expression itself (lines 275-276). What statistical test was used here? Did it consider that gene expression divergence is expected to be distributed between 0 and 1 while dominance divergence between -1 and 1? Please elaborate.

5. RNAseq data. Since poly-A selection was performed using total RNA, it should formally be mRNAseq data instead.

6. Figures are very informative. Fig. 1 & 7; I would suggest a more contrasting pattern of colors for increased clarity. Fig. 4 B; X-axis legend should be “Intersection class”, not “Interaction class”. Fig. 4C and D; the colors are helpful but I had difficulty seeing the asterisks of statistical significance so I would suggest increasing size a bit.

EDITORIAL COMMENTS

• Ln 103 “… single tissues, and/or highly diverged tissues….” replace with “… single tissues, and/or highly functionally diverged tissues….”

• Ln 132 “… From ASE data, we found that both cis and trans effects were more tissue-specific than genes with no differential expression…” replace with “… From ASE data, we found that genes with both cis and trans effects were more tissue-specific than genes with no differential expression… “

• Ln 151 “ … some overlap between SU58, ZI418 and their F1 hybrid and well as …” replace with “ … some overlap between SU58, ZI418 and their F1 hybrids, and as well as …”

• Ln 249 “In order to compare the magnitude of dominance regardless of which allele was dominant, we calculated h such that values between 0 and 1 or 0 and -1 represent varying degrees of additivity and dominance with values closer to 1 or -1 being more dominant and -1 representing complete dominance of the Swedish background and 1 representing complete dominance of the Zambian background, while values outside this range represent cases of overdominance of the respective background (see Methods for more details).” replace with “In order to compare the magnitude of dominance regardless of which allele was dominant, we calculated h such that values between 0 and 1 or 0 and -1 represent varying degrees of additivity and dominance, with values closer to -1 representing complete dominance of the Swedish background and 1 representing complete dominance of the Zambian background, while values outside this range represent cases of overdominance of the respective background (see Methods for more details).”

• Ln 351 “… and this was only significant after multiple test correction in the midgut…” replace with “… and this was significant after multiple test correction only in the midgut…”

• Ln 753 “ … the ASE analysis section above in order balance statistical power…” replace with “ … the ASE analysis section above in order to balance statistical power…”

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Genetics data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

10.1371/journal.pgen.1011257.r002
Author response to Decision Letter 0
Submission Version1
25 Jul 2024

Attachment Submitted filename: Response_to_Reviewers.pdf

10.1371/journal.pgen.1011257.r003
Decision Letter 1
Mackay Trudy F.C. Consulting Editor - PLoS Genetics
Williams Scott M. Section Editor
© 2024 Mackay, Williams
2024
Mackay, Williams
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version1
30 Jul 2024

Dear Dr. Glaser-Schmitt,

We are pleased to inform you that your manuscript entitled "Pervasive tissue-, genetic background-, and allele-specific gene expression effects in Drosophila melanogaster" has been editorially accepted for publication in PLOS Genetics. Congratulations!

Before your submission can be formally accepted and sent to production you will need to complete our formatting changes, which you will receive in a follow up email. Please be aware that it may take several days for you to receive this email; during this time no action is required by you. Please note: the accept date on your published article will reflect the date of this provisional acceptance, but your manuscript will not be scheduled for publication until the required changes have been made.

Once your paper is formally accepted, an uncorrected proof of your manuscript will be published online ahead of the final version, unless you’ve already opted out via the online submission form. If, for any reason, you do not want an earlier version of your manuscript published online or are unsure if you have already indicated as such, please let the journal staff know immediately at plosgenetics@plos.org.

In the meantime, please log into Editorial Manager at https://www.editorialmanager.com/pgenetics/, click the "Update My Information" link at the top of the page, and update your user information to ensure an efficient production and billing process. Note that PLOS requires an ORCID iD for all corresponding authors. Therefore, please ensure that you have an ORCID iD and that it is validated in Editorial Manager. To do this, go to ‘Update my Information’ (in the upper left-hand corner of the main menu), and click on the Fetch/Validate link next to the ORCID field.  This will take you to the ORCID site and allow you to create a new iD or authenticate a pre-existing iD in Editorial Manager.

If you have a press-related query, or would like to know about making your underlying data available (as you will be aware, this is required for publication), please see the end of this email. If your institution or institutions have a press office, please notify them about your upcoming article at this point, to enable them to help maximise its impact. Inform journal staff as soon as possible if you are preparing a press release for your article and need a publication date.

Thank you again for supporting open-access publishing; we are looking forward to publishing your work in PLOS Genetics!

Yours sincerely,

Trudy F.C. Mackay

Consulting Editor - PLoS Genetics

PLOS Genetics

Scott Williams

Section Editor

PLOS Genetics

www.plosgenetics.org

Twitter: @PLOSGenetics

----------------------------------------------------

Data Deposition

If you have submitted a Research Article or Front Matter that has associated data that are not suitable for deposition in a subject-specific public repository (such as GenBank or ArrayExpress), one way to make that data available is to deposit it in the Dryad Digital Repository. As you may recall, we ask all authors to agree to make data available; this is one way to achieve that. A full list of recommended repositories can be found on our website.

The following link will take you to the Dryad record for your article, so you won't have to re‐enter its bibliographic information, and can upload your files directly: 

http://datadryad.org/submit?journalID=pgenetics&manu=PGENETICS-D-24-00414R1

More information about depositing data in Dryad is available at http://www.datadryad.org/depositing. If you experience any difficulties in submitting your data, please contact help@datadryad.org for support.

Additionally, please be aware that our data availability policy requires that all numerical data underlying display items are included with the submission, and you will need to provide this before we can formally accept your manuscript, if not already present.

----------------------------------------------------

Press Queries

If you or your institution will be preparing press materials for this manuscript, or if you need to know your paper's publication date for media purposes, please inform the journal staff as soon as possible so that your submission can be scheduled accordingly. Your manuscript will remain under a strict press embargo until the publication date and time. This means an early version of your manuscript will not be published ahead of your final version. PLOS Genetics may also choose to issue a press release for your article. If there's anything the journal should know or you'd like more information, please get in touch via plosgenetics@plos.org.

10.1371/journal.pgen.1011257.r004
Acceptance letter
Mackay Trudy F.C. Consulting Editor - PLoS Genetics
Williams Scott M. Section Editor
© 2024 Mackay, Williams
2024
Mackay, Williams
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
17 Aug 2024

PGENETICS-D-24-00414R1

Pervasive tissue-, genetic background-, and allele-specific gene expression effects in Drosophila melanogaster

Dear Dr Glaser-Schmitt,

We are pleased to inform you that your manuscript entitled "Pervasive tissue-, genetic background-, and allele-specific gene expression effects in Drosophila melanogaster" has been formally accepted for publication in PLOS Genetics! Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out or your manuscript is a front-matter piece, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Genetics and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Zsofia Freund

PLOS Genetics

On behalf of:

The PLOS Genetics Team

Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom

plosgenetics@plos.org | +44 (0) 1223-442823

plosgenetics.org | Twitter: @PLOSGenetics
==== Refs
References

1 Brawand D , Soumillon M , Necsulea A , Julien P , Csárdi G , Harrigan P , et al . The evolution of gene expression levels in mammalian organs. Nature. 2011;478 (7369 ):343–8. doi: 10.1038/nature10532 22012392
2 Consortium G. Human genomics. The Genotype-Tissue Expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science. 2015;348 (6235 ):648–60. doi: 10.1126/science.1262110 25954001
3 Shalek AK , Satija R , Shuga J , Trombetta JJ , Gennert D , Lu D , et al . Single-cell RNA-seq reveals dynamic paracrine control of cellular variation. Nature. 2014;510 (7505 ):363–9. doi: 10.1038/nature13437 24919153
4 Witt E , Benjamin S , Svetec N , Zhao L . Testis single-cell RNA-seq reveals the dynamics of de novo gene transcription and germline mutational bias in Drosophila. Elife. 2019;8.
5 Coolon JD , McManus CJ , Stevenson KR , Graveley BR , Wittkopp PJ . Tempo and mode of regulatory evolution in Drosophila. Genome Res. 2014;24 (5 ):797–808.24567308
6 Chen J , Nolte V , Schlötterer C . Temperature stress mediates decanalization and dominance of gene expression in Drosophila melanogaster. PLoS Genet. 2015;11 (2 ):e1004883.25719753
7 King MC , Wilson AC . Evolution at two levels in humans and chimpanzees. Science. 1975;188 (4184 ):107–16. doi: 10.1126/science.1090005 1090005
8 Wray GA , Hahn MW , Abouheif E , Balhoff JP , Pizer M , Rockman MV , et al . The Evolution of Transcriptional Regulation in Eukaryotes. Molecular Biology and Evolution. 2003;20 (9 ):1377–419. doi: 10.1093/molbev/msg140 12777501
9 Buchberger E , Reis M , Lu TH , Posnien N . Cloudy with a Chance of Insights: Context Dependent Gene Regulation and Implications for Evolutionary Studies. Genes (Basel). 2019;10 (7 ). doi: 10.3390/genes10070492 31261769
10 Ayroles JF , Carbone MA , Stone EA , Jordan KW , Lyman RF , Magwire MM , et al . Systems genetics of complex traits in Drosophila melanogaster. Nat Genet. 2009;41 (3 ):299–307.19234471
11 Mackay TF , Stone EA , Ayroles JF . The genetics of quantitative traits: challenges and prospects. Nat Rev Genet. 2009;10 (8 ):565–77. doi: 10.1038/nrg2612 19584810
12 Barbeira AN , Dickinson SP , Bonazzola R , Zheng J , Wheeler HE , Torres JM , et al . Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics. Nat Commun. 2018;9 (1 ):1825. doi: 10.1038/s41467-018-03621-1 29739930
13 Signor SA , Nuzhdin SV . The Evolution of Gene Expression in cis and trans. Trends Genet. 2018;34 (7 ):532–44.29680748
14 Hill MS , Vande Zande P , Wittkopp PJ . Molecular and evolutionary processes generating variation in gene expression. Nat Rev Genet. 2021;22 (4 ):203–15. doi: 10.1038/s41576-020-00304-w 33268840
15 Wittkopp PJ , Haerum BK , Clark AG . Evolutionary changes in cis and trans gene regulation. Nature. 2004;430 (6995 ):85–8.15229602
16 Wittkopp PJ , Haerum BK , Clark AG . Regulatory changes underlying expression differences within and between Drosophila species. Nat Genet. 2008;40 (3 ):346–50.18278046
17 Graze RM , McIntyre LM , Main BJ , Wayne ML , Nuzhdin SV . Regulatory divergence in Drosophila melanogaster and D. simulans, a genomewide analysis of allele-specific expression. Genetics. 2009;183 (2 ):547–61, 1SI-21SI.19667135
18 Metzger BPH , Wittkopp PJ , Coolon JD . Evolutionary Dynamics of Regulatory Changes Underlying Gene Expression Divergence among Saccharomyces Species. Genome Biol Evol. 2017;9 (4 ):843–54.28338820
19 Glaser-Schmitt A , Zečić A , Parsch J . Gene Regulatory Variation in Drosophila melanogaster Renal Tissue. Genetics. 2018;210 (1 ):287–301.29976765
20 McManus CJ , Coolon JD , Duff MO , Eipper-Mains J , Graveley BR , Wittkopp PJ . Regulatory divergence in Drosophila revealed by mRNA-seq. Genome Res. 2010;20 (6 ):816–25.20354124
21 Osada N , Miyagi R , Takahashi A . Cis—and Trans -regulatory Effects on Gene Expression in a Natural Population of Drosophila melanogaster. Genetics. 2017;206 (4 ):2139–48.28615283
22 Benowitz KM , Coleman JM , Allan CW , Matzkin LM . Contributions of cis- and trans-Regulatory Evolution to Transcriptomic Divergence across Populations in the Drosophila mojavensis Larval Brain. Genome Biol Evol. 2020;12 (8 ):1407–18.32653899
23 Majane AC , Cridland JM , Blair LK , Begun DJ . Evolution and genetics of accessory gland transcriptome divergence between Drosophila melanogaster and D. simulans. Genetics. 2024.
24 Verta JP , Jones FC . Predominance of cis-regulatory changes in parallel expression divergence of sticklebacks. Elife. 2019;8. doi: 10.7554/eLife.43785 31090544
25 Durkin SM , Ballinger MA , Nachman MW . Tissue-specific and cis-regulatory changes underlie parallel, adaptive gene expression evolution in house mice. PLoS Genet. 2024;20 (2 ):e1010892.38306396
26 Fear JM , León-Novelo LG , Morse AM , Gerken AR , Van Lehmann K , Tower J , et al . Buffering of Genetic Regulatory Networks in Drosophila melanogaster. Genetics. 2016;203 (3 ):1177–90.27194752
27 Meiklejohn CD , Coolon JD , Hartl DL , Wittkopp PJ . The roles of cis- and trans-regulation in the evolution of regulatory incompatibilities and sexually dimorphic gene expression. Genome Res. 2014;24 (1 ):84–95.24043293
28 Puixeu G , Macon A , Vicoso B . Sex-specific estimation of cis and trans regulation of gene expression in heads and gonads of Drosophila melanogaster. G3 (Bethesda). 2023;13 (8 ).
29 Kapopoulou A , Kapun M , Pieper B , Pavlidis P , Wilches R , Duchen P , et al . Demographic analyses of a new sample of haploid genomes from a Swedish population of Drosophila melanogaster. Sci Rep. 2020;10 (1 ):22415. doi: 10.1038/s41598-020-79720-1 33376238
30 Pool JE , Corbett-Detig RB , Sugino RP , Stevens KA , Cardeno CM , Crepeau MW , et al . Population Genomics of sub-saharan Drosophila melanogaster: African diversity and non-African admixture. PLoS Genet. 2012;8 (12 ):e1003080.23284287
31 Sprengelmeyer QD , Mansourian S , Lange JD , Matute DR , Cooper BS , Jirle EV , et al . Recurrent Collection of Drosophila melanogaster from Wild African Environments and Genomic Insights into Species History. Mol Biol Evol. 2020;37 (3 ):627–38.31730190
32 González J , Macpherson JM , Petrov DA . A recent adaptive transposable element insertion near highly conserved developmental loci in Drosophila melanogaster. Mol Biol Evol. 2009;26 (9 ):1949–61.19458110
33 Mateo L , Ullastres A , González J . A transposable element insertion confers xenobiotic resistance in Drosophila. PLoS Genet. 2014;10 (8 ):e1004560. doi: 10.1371/journal.pgen.1004560 25122208
34 Glaser-Schmitt A , Parsch J . Functional characterization of adaptive variation within a cis-regulatory element influencing Drosophila melanogaster growth. PLoS Biol. 2018;16 (1 ):e2004538.29324742
35 Ramnarine TJS , Glaser-Schmitt A , Catalán A , Parsch J . Population Genetic and Functional Analysis of a cis-Regulatory Polymorphism in the Drosophila melanogaster Metallothionein A gene. Genes. 2019;10 (2 ):147.30769915
36 Glaser-Schmitt A , Parsch J . Dynamics and stage-specificity of between-population gene expression divergence in the Drosophila melanogaster larval fat body. PLoS Genet. 2023;19 (4 ):e1010730.37099612
37 Miguel-Aliaga I , Jasper H , Lemaitre B . Anatomy and Physiology of the Digestive Tract of Drosophila melanogaster. Genetics. 2018;210 (2 ):357–96.30287514
38 Cohen E , Sawyer JK , Peterson NG , Dow JAT , Fox DT . Physiology, Development, and Disease Modeling in the Drosophila Excretory System. Genetics. 2020;214 (2 ):235–64. doi: 10.1534/genetics.119.302289 32029579
39 Liu Q , Jin LH . Organ-to-Organ Communication: A Drosophila Gastrointestinal Tract Perspective. Front Cell Dev Biol. 2017;5 :29. doi: 10.3389/fcell.2017.00029 28421183
40 Lemaitre B , Miguel-Aliaga I . The digestive tract of Drosophila melanogaster. Annu Rev Genet. 2013;47 :377–404.24016187
41 Krause SA , Overend G , Dow JAT , Leader DP . FlyAtlas 2 in 2022: enhancements to the Drosophila melanogaster expression atlas. Nucleic Acids Res. 2022;50 (D1 ):D1010–D5.34718735
42 Lavington E , Kern AD . The Effect of Common Inversion Polymorphisms In(2L) and In(3R)Mo on Patterns of Transcriptional Variation in Drosophila melanogaster. G3 (Bethesda). 2017;7 (11 ):3659–68.28916647
43 Corbett-Detig RB , Hartl DL . Population genomics of inversion polymorphisms in Drosophila melanogaster. PLoS Genet. 2012;8 (12 ):e1003056.23284285
44 Lack JB , Lange JD , Tang AD , Corbett-Detig RB , Pool JE . A Thousand Fly Genomes: An Expanded Drosophila Genome Nexus. Mol Biol Evol. 2016;33 (12 ):3308–13. doi: 10.1093/molbev/msw195 27687565
45 Said I , Byrne A , Serrano V , Cardeno C , Vollmers C , Corbett-Detig R . Linked genetic variation and not genome structure causes widespread differential expression associated with chromosomal inversions. Proc Natl Acad Sci U S A. 2018;115 (21 ):5492–7. doi: 10.1073/pnas.1721275115 29735663
46 Broderick NA , Buchon N , Lemaitre B . Microbiota-induced changes in Drosophila melanogaster host gene expression and gut morphology. mBio. 2014;5 (3 ):e01117–14.
47 Wong CN , Ng P , Douglas AE . Low-diversity bacterial community in the gut of the fruitfly Drosophila melanogaster. Environ Microbiol. 2011;13 (7 ):1889–900.21631690
48 Chandler JA , Lang JM , Bhatnagar S , Eisen JA , Kopp A . Bacterial communities of diverse Drosophila species: ecological context of a host-microbe model system. PLoS Genet. 2011;7 (9 ):e1002272. doi: 10.1371/journal.pgen.1002272 21966276
49 Hung RJ , Hu Y , Kirchner R , Liu Y , Xu C , Comjean A , et al . A cell atlas of the adult Drosophila midgut. Proc Natl Acad Sci U S A. 2020;117 (3 ):1514–23.31915294
50 Costa M , Reeve S , Grumbling G , Osumi-Sutherland D . The Drosophila anatomy ontology. J Biomed Semantics. 2013;4 (1 ):32. doi: 10.1186/2041-1480-4-32 24139062
51 Carroll SB . Endless Forms: The Evolution of Gene Regulation and Morphological Diversity. Cell. 2000;101 (6 ):577–80. doi: 10.1016/s0092-8674(00)80868-5 10892643
52 Carroll SB . Evo-Devo and an Expanding Evolutionary Synthesis: A Genetic Theory of Morphological Evolution. Cell. 2008;134 (1 ):25–36. doi: 10.1016/j.cell.2008.06.030 18614008
53 Stern DL , Orgogozo V . The loci of evolution: how predictable is genetic evolution? Evolution. 2008;62 (9 ):2155–77. doi: 10.1111/j.1558-5646.2008.00450.x 18616572
54 Catalán A , Glaser-Schmitt A , Argyridou E , Duchen P , Parsch J . An Indel Polymorphism in the MtnA 3’ Untranslated Region Is Associated with Gene Expression Variation and Local Adaptation in Drosophila melanogaster. PLoS Genet. 2016;12 (4 ):e1005987.27120580
55 Daborn PJ , Yen JL , Bogwitz MR , Le Goff G , Feil E , Jeffers S , et al . A Single P450 Allele Associated with Insecticide Resistance in Drosophila. Science. 2002;297 (5590 ):2253–6.12351787
56 Simhadri RK , Fast EM , Guo R , Schultz MJ , Vaisman N , Ortiz L , et al . The Gut Commensal Microbiome of Drosophila melanogaster Is Modified by the Endosymbiont Wolbachia. mSphere. 2017;2 (5 ). doi: 10.1128/mSphere.00287-17 28932814
57 Staubach F , Baines JF , Künzel S , Bik EM , Petrov DA . Host species and environmental effects on bacterial communities associated with Drosophila in the laboratory and in the natural environment. PLoS One. 2013;8 (8 ):e70749. doi: 10.1371/journal.pone.0070749 23967097
58 Wang Y , Kapun M , Waidele L , Kuenzel S , Bergland AO , Staubach F . Common structuring principles of the Drosophila melanogaster microbiome on a continental scale and between host and substrate. Environ Microbiol Rep. 2020;12 (2 ):220–8.32003146
59 Coolon JD , Stevenson KR , McManus CJ , Graveley BR , Wittkopp PJ . Genomic imprinting absent in Drosophila melanogaster adult females. Cell Rep. 2012;2 (1 ):69–75.22840398
60 Callahan BJ , McMurdie PJ , Rosen MJ , Han AW , Johnson AJ , Holmes SP . DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13 (7 ):581–3. doi: 10.1038/nmeth.3869 27214047
61 Yilmaz P , Parfrey LW , Yarza P , Gerken J , Pruesse E , Quast C , et al . The SILVA and "All-species Living Tree Project (LTP)" taxonomic frameworks. Nucleic Acids Res. 2014;42 (Database issue ):D643–8. doi: 10.1093/nar/gkt1209 24293649
62 Wickham H. ggplot2: Elegant Graphics for Data Analysis: Springer-Verlag New York; 2016.
63 McMurdie PJ , Holmes S . phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One. 2013;8 (4 ):e61217. doi: 10.1371/journal.pone.0061217 23630581
64 Oksanen J , Simpson GL , Blanchet FG , Kindt R , Legendre P , Minchin PR , et al . vegan: Community Ecology Package. 2022.
65 Bates D , Mächler M , Bolker B , Walker S . Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software. 2015;67 (1 ):1–48.
66 Lenth RV , Bolker B , Buerkner P , Giné-Vázquez I , Herve M , Jung M , et al . emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.10.0 ed 2024.
67 Lack JB , Cardeno CM , Crepeau MW , Taylor W , Corbett-Detig RB , Stevens KA , et al . The Drosophila genome nexus: a population genomic resource of 623 Drosophila melanogaster genomes, including 197 from a single ancestral range population. Genetics. 2015;199 (4 ):1229–41.25631317
68 Thurmond J , Goodman JL , Strelets VB , Attrill H , Gramates LS , Marygold SJ , et al . FlyBase 2.0: the next generation. Nucleic Acids Res. 2019;47 (D1 ):D759–D65. doi: 10.1093/nar/gky1003 30364959
69 Sedlazeck FJ , Rescheneder P , von Haeseler A . NextGenMap: fast and accurate read mapping in highly polymorphic genomes. Bioinformatics. 2013;29 (21 ):2790–1. doi: 10.1093/bioinformatics/btt468 23975764
70 Love MI , Huber W , Anders S . Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15 (12 ):550. doi: 10.1186/s13059-014-0550-8 25516281
71 Yanai I , Benjamin H , Shmoish M , Chalifa-Caspi V , Shklar M , Ophir R , et al . Genome-wide midrange transcription profiles reveal expression level relationships in human tissue specification. Bioinformatics. 2005;21 (5 ):650–9. doi: 10.1093/bioinformatics/bti042 15388519
72 Falconer DS , Mackay TF . Introduction to Quantitative Genetics. 4th ed. Harlow: Pearson Education; 1995.
73 Glaser-Schmitt A , Wittmann MJ , Ramnarine TJS , Parsch J . Sexual Antagonism, Temporally Fluctuating Selection, and Variable Dominance Affect a Regulatory Polymorphism in Drosophila melanogaster. Mol Biol Evol. 2021;38 (11 ):4891–907.34289067
74 Gibson G , Riley-Berger R , Harshman L , Kopp A , Vacha S , Nuzhdin S , et al . Extensive sex-specific nonadditivity of gene expression in Drosophila melanogaster. Genetics. 2004;167 (4 ):1791–9.15342517
75 Benjamini Y , Hochberg Y . Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B. 1995;57 (1 ):289–300.
76 Smith RN , Aleksic J , Butano D , Carr A , Contrino S , Hu F , et al . InterMine: a flexible data warehouse system for the integration and analysis of heterogeneous biological data. Bioinformatics. 2012;28 (23 ):3163–5. doi: 10.1093/bioinformatics/bts577 23023984
