Contrasting patterns of divergence at the regulatory and sequence level in European Daphnia galeata natural populations

Abstract

Understanding the genetic basis of local adaptation has long been a focus of evolutionary biology. Recently, there has been increased interest in deciphering the evolutionary role of Daphnia's plasticity and the molecular mechanisms of local adaptation. Using transcriptome data, we assessed the differences in gene expression profiles and sequences in four European Daphnia galeata populations. In total, ~33% of 32,903 transcripts were differentially expressed between populations. Among 10,280 differentially expressed transcripts, 5,209 transcripts deviated from neutral expectations and their population-specific expression pattern is likely the result of local adaptation processes. Furthermore, a SNP analysis allowed inferring population structure and distribution of genetic variation. The population divergence at the sequence level was comparatively higher than the gene expression level by several orders of magnitude consistent with strong founder effects and lack of gene flow between populations. Using sequence homology, the candidate transcripts were annotated using a comparative genomics approach. Additionally, we also performed a weighted gene co-expression analysis to identify population-specific regulatory patterns of transcripts in D. galeata. Thus, we identified candidate transcriptomic regions for local adaptation in this key species of aquatic ecosystems in the absence of any laboratory-induced stressor.

Keywords: DRIFTSEL; RNA‐seq; WGCNA; constitutive gene expression; molecular phenotype; population transcriptomics.

Figure 1

Waterfleas, Daphnia galeata, Photograph: Mathilde Cordellier

Figure 2

Gene expression patterns. (a) Gene expression PCA of the four sampled populations: Pop.G (Lake Greifensee), Pop.J (Jordan Reservoir), Pop.LC (Lake Constance), and Pop.M (Müggelsee). Percentages on the x‐ and y‐axis indicate the percentage of variance explained by each principal component. (b) Venn diagram illustrating the number of differentially expressed transcripts (DET) between the four populations. Numbers in brackets indicate the number of transcripts deviating from the neutral expectations according to the DRIFTSEL analysis

Figure 3

SNP patterns and inbreeding coefficient. (a) SNP PCA of the four sampled populations: Pop.G (Lake Greifensee), Pop.J (Jordan reservoir), Pop.LC (Lake Constance), and Pop.M (Müggelsee). Percentages on the x‐ and y‐axis indicate the percentage of variance explained by each principal component. (b) Barplot illustrating the inbreeding coefficient for each genotype

Figure 4

Differentiating misassembly from inparalogs. (a) Barplot showing the number of DETs co‐occurring with DETs from other populations within an orthoMCL cluster. “No‐cluster DETs” refers to DETs not assigned to an orthoMCL cluster. 1Pop, 2Pop, 3Pop, and 4Pop refer to DETs found in orthoMCL clusters containing at least one, two, three, and four population(s), respectively. (b) Histogram of pairwise sequence divergence values calculated for all D. galeata sequences co‐occurring in an orthoMCL cluster belonging to 2Pop, 3Pop, and 4Pop categories

Figure 5

Flow diagram representing the proportion of transcripts that are candidates for local adaptation at the regulatory and sequence level. Each analysis or “step” is represented by a vertical group of black rectangle bars, called nodes. The colored areas linking the nodes are called “flows.” The DESeq2 step contains four nodes: PopG (yellow), PopJ (black), PopLC (pink), and PopM (green), which represent the number of transcripts specifically upregulated in each of the four populations as identified by DESeq2 analysis. The DRIFTSEL step contains 2 nodes: “H.value ≤0.95” (grey) and “H.value ≥0.95” (purple). The LOSITAN step contains 5 nodes: “NC” (grey) with transcripts without LOSITAN result (not calculated); “noOL” (grey), transcripts where none of the SNPs in a transcript were identified as outliers; “Bal” (cyan), transcripts containing at least one SNP that is under balancing selection; “Div” (pink), transcripts containing at least one SNP under diversifying selection; and “BalDiv” (pale green), transcripts containing SNPs that are under both balancing and diversifying selection. The Tajima’s D step contains eight nodes. Each node classifies the transcripts according to the obtained Tajima’s D values. “AllNeg” means that transcripts have a negative D value in all four populations; “AllPos” means that transcripts have a positive D value in all four populations; “AllNonSig” means transcripts have nonsignificant D values in all four populations; “NegNonsig” means transcripts in the four populations have either a negative D value or a nonsignificant D value; “PosNonsig” means transcripts in the four populations have either a positive D value or a nonsignificant D value; “PosNeg” means transcripts in the four populations have either a positive or negative D value; “PosNegNonsig” means transcripts in the four populations have either a positive or negative or an insignificant D value

Figure 6

Cluster dendrogram of transcripts for the reference network in Daphnia galeata, with dissimilarity based on the topological overlap matrix (TOM). The co‐expression modules are colored in an arbitrary way by the WGCNA package, and the size of the bar is proportional to the number of transcripts in the module. The right‐hand side grid represents the module conservation in each population. Modules with a Z‐score ≤10 are shown in white, and modules with a Z‐score ≥10 are colored in dark grey