Gene-environment and protein-degradation signatures characterize genomic and phenotypic diversity in wild Caenorhabditis eleganspopulations

Local C. eleganspopulations are genotypically separable

Previous investigations have studied population genetics and genomic diversity in C. elegans, focusing on global [31–35] or local [17, 36] populations. Cutter showed that there is a lack of geographic distribution of C. elegans genome sequences [31], and Andersen et al. reported that chromosome-scale selective sweeps have acted to reduce genetic variation, and have shaped the global C. elegans population structure in recent history [37]. Barrière and Félix concluded that local diversity in this organism is high [17]. In all of these papers, diversity in C. elegans was measured as genetic diversity. To date, very few papers have been published concerning phenotypic variation in wild isolates (by this we mean isolates that have not been maintained in the laboratory for a long time), and studies used only a small number of isolates [38–40]. To our knowledge, no studies have been reported on genotype-phenotype relations in wild populations.

As mentioned above, we used wild strains from Santeuil and Orsay and a number of out-group strains. All these were genotyped based on the hybridization of genomic C. elegans DNA to microarrays (see Methods section for details), resulting in the identification of 6,368 polymorphic genes with an absolute ratio of 0.5 for the mean hybridization intensity (Figure 1A). Most of these (around 66%) were found in only one to three strains (Figure 1B), showing that between-strain variation is more abundant than between-site variation. Major hotspots of polymorphic genes were found on chromosomes II and V and minor on the other chromosomes (Figure 1C; see Additional file 1 worksheet B; see Additional files 2, 3, 4). The hotspots of polymorphic genes co-localize with the c-type lectin, nuclear hormone receptor (nhr), and math genes, and chemoreceptor gene clusters [41]. In addition, microsatellite loci were used to determine the population structure [36] of the Orsay and Santeuil strains (see Additional file 1, worksheet C; see Additional file 5, panel A).

The two isolation sites were genotypically separable. This was shown by analyses of the presence of gene polymorphisms using principal component analysis (PCA) (Figure 2A), a distance matrix visualized by an unrooted neighbor-joining (NJ) tree (Figure 2B), and a minimum spanning network of the microsatellite data (see Additional file 5, panel A). The minimum spanning network, PCA, and NJ analyses showed a clear distinction between the Santeuil and the Orsay strains, with one large genetic group and several smaller genetic groups being identified for both isolation sites. In PCA, the first two principal components capture around 75% of the variation in DNA hybridizations. As shown in Figure 2A, N2 is in the far right top corner, indicating its genetic difference from all other strains. Moreover, the NJ tree showed that the Santeuil strains (groups S1, to S3) and Orsay strains are different from both N2 and CB4856. In both the PCA and NJ analyses, the Orsay group (group O) was seen to be genetically less diverse than the Santeuil group. Within the main Santeuil group (group S1: all Santeuil strains except JU1924, JU1925, JU1926, JU1934, JU1935, and JU1936), diversity was slightly larger. Furthermore, in the NJ tree two small genotypic groups were separate from the main Orsay and Santeuil groups (group S2: JU1924, JU1925, and JU1926; group S3: JU1511, JU1934, JU1935, and JU1936). The strains within these separate groups were all from the Santeuil site, except for JU1511, which is from the Orsay site. The strains from Santeuil in group S2 were isolated from a single hogweed stem. Similarly, those from S3 were also isolated from their own single hogweed stem (see Additional file 1, worksheet A). Other strains were found on different hogweed stems. We found that strains isolated from an individual hogweed stem grouped close together, but were not found to form their own separate genotypic groups.

The Dutch strains were isolated from two isolation sites, with WN2001 isolated from one site and WN2002 and WN2003 from the other. The latter two strains grouped together in the NJ tree. One of the French out-group strains, JU396, grouped with the Santeuil strains, whereas the other, JU314, was different from the rest of the strains. N2 and CB4856 were as diverse from the other out-group members as the whole out-group was from the Orsay or Santeuil strains. By comparing the four genetic groups (O, S1, S2, and S3) with the out-group, the genes that were polymorphic were identified (see Additional file 5, panel B). Group S1 appeared to be the most divergent from the out-group, with 3,181 genes that differed significantly (false discovery rate (FDR) = 0.05).

The genetic separation between the Orsay and Santeuil populations was prevalent on all chromosomes (Figure 3). From the scale of the axes, it can be seen that most chromosomes contributed to the separation between the two isolation sites and the out-group, except for chromosome II; on this chromosome the Santeuil and Orsay lines formed one group that was separate only from the out-group. Chromosome V was the largest contributor to the variation between Orsay and Santeuil, most likely because of the generally higher level of variation among these strains (see Additional file 4). Of the approximately 2,500 genes that were different between S1, S2, and S3, around 1,050 are located on chromosome V. Of all the genes on chromosome V, around 210 are polymorphic between S1, S2, and S3. This is a significant enrichment (P<1 × 10^-76) when compared with the other chromosomes, of which 8 to 10% of the genes are polymorphic. Chromosomes I and X are under-represented for polymorphic genes (both around 8% and P<1 × 10^-12).

The detected genotypic diversity between sites is in line with genotyping results from 31 markers using amplified fragment length polymorphism analysis, and with microsatellite results from two loci in other local C. elegans populations [17]. However, selective sequencing using restriction site-associated DNA tags did not reveal significant local diversity, possibly because of the limited number of different genotypes per location [37].

The four genetic groups (O, S1, S2, and S3) identified by analysis of the ratio intensities were used as input to search for all the genes linked to each isolation site. This allowed us to identify polymorphic genes by minor hybridization differences (absolute ratio (AR) <0.5), beyond those 6,386 found by major hybridization differences (AR >0.5). In this way, we identified 3,742 genes (FDR = 0.05) that were linked to each isolation site (see Additional file 6, panel A). Of these 3,742 genes, 2,403 were already identified as highly polymorphic in the initial analysis on major hybridization differences, and an additional 1,339 genes associated with isolation were found with only minor hybridization differences. Of the genes with major hybridization differences, around 62% could not be linked to an isolation site, again showing that between-strain variation is more abundant than between-site variation. Genes linked to isolation site were fond to be enriched for the gene classes fbox, math, bath, btb, C-type lectin (clec), serpentine chemoreceptor, and nhr.

Local C. eleganspopulations are separable on the basis of their transcriptomes

Next, the influence of natural genetic variation on gene expression was studied by measuring the transcript levels of all genes of all strains, corrected for differential hybridization. An NJ tree was constructed based on the RNA hybridization data (Figure 2C). This tree showed that the genetic groups O, S1, S2, S3, and the out-group, were also separable based on gene expression level. Again, CB4856 and N2 differed from most of the other natural strains. Isolation site and genetic group influenced the variation in RNA levels of 6,930 and 7,996 genes, respectively (see Additional file 7). Most of these genes (77% and 78%) were not influenced by DNA polymorphisms (which were the cause of variation in 2,330 genes) or genotype (affecting 773 or 1,336 genes, depending on the genetic group or isolation site that was incorporated in the model).

Expressed genes linked to isolation site are enriched for the gene classes clec, fbxa, bath, and nhr

Enrichment analyses were performed for the genes whose RNA levels were influenced by isolation site or genetic group (see Additional file 1, worksheets L–N). The gene classes clec, fbxa, bath, and nhr were significantly enriched, thus yielding similar results to our DNA-level enrichment analyses. In addition, several nematode-specific peptide families were also enriched.

Together, these results show that at the genomic level, variation between local populations is enriched for the gene classes fbox, math, bath, btb, clec, serpentine chemoreceptor, and nhr, many of which are involved in gene-environment interactions [42–45]. Interestingly, we found that the gene classes clec, fbxa, bath, and nhr were also enriched, with variation linked to isolation site on the transcriptional level, even though the strains that originated from different sites were cultured under the same conditions. It has been shown that many of these groups of genes are differentially expressed after pathogen exposure, and thus could be involved in the immune response. For instance, C-type lectin domain-containing proteins (CTLD proteins, gene class clec) have been repeatedly proposed to contribute to nematode immunity [46]. The immune function of these genes is supported by their specific upregulation in infected C. elegans[43, 44, 47–50] and also by reduced immune phenotypes after RNAi knockdown of clec-70, clec-17, clec-60, or clec-86[7, 51]. Furthermore, F-box proteins (gene class fbxa) are part of the protein degradation pathway [52]. In this pathway, substrates for degradation are ubiquinated to be recognizable by the 26S proteasome. Taken together, these results show that local genetic diversity reflects site-specific signatures of immune response and protein degradation pathways in C. elegans. We also found that, in addition to genotypes, transcript profiles can be used to distinguish between local C. elegans populations, and may indicate the functional importance of the identified genes or gene classes in different environments [53].

Polymorphic genes are enriched for genes involved in gene-environment interactions

Polymorphic genes between the populations were compared using a set of genes extracted from 19 different studies on gene expression in C. elegans exposed to biotic and abiotic factors (see Additional file 1, worksheet S). In the wild, C. elegans is exposed to many different bacteria. In studies on the effect on gene expression of various bacteria, such as Lactobacillus rhamnosus[54], Microbacterium. nematophilum[7], Drechmeria coniospora[55], Serratia marcescens[8], Xenorhabdus nematophila[8], and Pseudomonas aeruginosa[56], c-type lectins were always found to be differentially expressed, as were in most cases the F-box protein genes. Receptors that are used to sense the environment, such as nhrs and serpentine receptors are also frequently differentially expressed when C. elegans is exposed to different bacteria. In response to abiotic factors such as temperature [12], osmotic stress [10] or ions [57, 58], the c-type lectins and F-box protein genes are also always differentially expressed. Furthermore, the c-type lectins, F-box protein genes, and receptor genes are differentially expressed in the presence of various other substances that can be encountered by wild C. elegans strains, including tryptophan [59], β-naphthoflavone [60], H₂S [61], fluoranthene [62], hormones [63], sediment [64], humic substances [65], and pesticides [12, 66, 67]. The other gene classes (bath, math, and btb) that are importaty for the variation seen between the locations at which the wild C. elegans strains were isolated were also found to be differentially expressed in several of the aforementioned environmental studies. Altogether, the differential expression of genes in environmental studies indicates that the genes that are important for the variation between local populations of C. elegans are indeed of significance for interactions with the environment.

Local populations are separable for some fitness traits

Local populations are separable with regard to food preference

We then investigated if the wild strains differed from each other in their food preference behavior for naturally co-existing bacteria and for E. coli. Under the laboratory conditions we used, the worms preferred E. coli OP50 over all other bacteria, followed by Erwinia rhapontici, Sphingobacterium sp., Rhodococcus erythropolis, and Lactococcus lactis (Figure 4; see Additional file 1, worksheet P). Worms from Santeuil preferred E. rhapontici (isolated from Santeuil) equally to E. coli., whereas worms from Orsay preferred E. coli over E. rhapontici. This suggests that Santeuil worms could have a slight preference for the bacterium species with which they are more likely to be familiar for an overview of the average preference of all strains, see Additional file 9B; for the complete dataset, see Additional file 1, worksheet Q (Wormcount and Choice Index) and see Additional file 1, worksheet R (Significances). The most significant differences found between the Orsay and Santeuil strains were for the bacterial combinations E. coli OP50/E. rhapontici, E. coli OP50/R. erythropolis. and E. rhapontici/Sphingobacterium sp. (Figure 4).

The canonical strains Bristol N2 and CB4856 are genetic outliers

We finally investigated how far the canonical wild types Bristol N2 and CB4856 relate to the recently isolated natural strains. Both genotypic and transcriptomic analyses identified Bristol N2 and CB4856 as clearly differing from the standing genetic variation of the wild isolated strains. This also applied when we compared N2 and CB4856 with some of the other out-group strains. It is assumed that many alleles in CB4856 and N2 are laboratory-derived because both N2 and CB4856 went through multiple phenotypic and genetic bottlenecks over the past decades of laboratory maintenance [1, 68]. Together with our results, this shows that the genotype-phenotypes relations in both N2 and CB4856 are likely to be obscured by a number of laboratory-derived alleles with large effect. This might impede the detection and functional analysis of many other genes and alleles that by themselves have small effects on phenotypes, but together might have a large effect.

Nematode and bacterial strains

The main set of strains of C. elegans comprised 41 new strains that were isolated (by M-A Félix) from two different locations in France (Orsay and Santeuil). The out-group comprised three new strains isolated in the Netherlands, two strains previously isolated in France, and the most diverse canonical strains N2 (Bristol) and CB4856 (Hawaii) [16, 18–21, 23–29, 67] (see Additional file 1, worksheet A for details). All strains were routinely maintained on nematode growth medium (NGM) with E. coli OP50 as a food source [69]. E. coli OP50 was used in all experiments, except for the population growth experiment, in which B. thuringiensis NRRL B-18247 and B. thuringiensis DSM-350 were used next to E. coli[70]. In the food preference experiment, in addition to E. coli OP50, E. rhapontici and R.s erythropolis (both isolated from and unique for Santeuil), and L. lactis and Sphingobacterium sp. (both isolated from and unique for Orsay) were used (all bacteria were isolated and identified by M-A. Félix and B. Samuel).

Genomic DNA analysis: worm culturing, DNA isolation, DNA microarrays, and statistical analysis

where I is the DNA hybridization intensity of genotype x, G is the genetic group (out-group, S1, S2, or S3) and E is the error. The thresholds used were –log10(p) 2.5, 2.3, 3.2, and 3.3 for O, S1, S2, and S3 respectively (for the number of genes per genetic group and overlapping genes (see Additional file 5B). The significance thresholds, adjusted for multiple testing, were determined by permutation, for which the same model was used, with the spot intensities randomly distributed over the genotypes (a P value that gave a ratio of false positives/true positives of <0.05 was used).

mRNA analysis: culturing, isolation, RNA-microarrays, and statistical analysis

For the mRNA microarrays, any males were discarded and only hermaphrodites grown on E. coli OP50 were used. Two independent replicates of each strain (synchronized late L4 larvae) were analyzed. For mRNA isolation, a commercial kit (RNEasy Micro Kit; Qiagen, Hilden, Germany) was used, following the manufacturer’s protocol (Purification of Total RNA from Animal and Human Tissues) provided with the kit, with modified lysing procedure (see Additional file 10). The microarrays used were as described above (C. elegans (V2) Gene Expression Microarray 4X44K slides; Agilent) following the manufacturer’s instructions. For processing of the RNA microarray data, the limma package for the R software environment was used. No background correction was performed, as recommended previously [71]. For within-array normalization, the LOESS method was used and for between-array normalization, the quantile method was used. Expression variation was determined by linear models. The variation in intensities could be explained by batch, DNA hybridization, genetic group, and genotype (see also the paragraph on statistics in the Genomic DNA Analysis section). Significance thresholds, adjusted for multiple testing, were determined by permutations of all spots on the array. In the permutations, the RNA hybridization intensities were randomly distributed over the genotypes and batches (the P-value that gave a ratio of false positives/true positives of <0.05 was used).

Enrichment analysis

All enrichment analyses were performed using a hypergeometric test. The number of genes selected by a criterion in this paper (for example, linked to a genetic group) were compared with the genes with a specific annotation (for example, c-type lectin). The chance that a number of genes will be overlapping depends on the total group size, the number of genes selected, and the number of genes with a specific annotation. This chance, together with the number of overlapping genes, can be used in a hypergeometric test. Annotation groups were considered enriched when the overlap was more than three genes and the significance –log10(p) was greater than 2.5.

Polymorphic genes between populations were compared with a set of differentially expressed genes extracted from a diverse set of gene-environment interaction studies in C. elegans. All enrichment analyses were performed using a hypergeometric test.

Phenotypic assays

Microsatellite analysis

Population genetic differentiation was assessed using six microsatellite loci (see Additional file 1, worksheet C), which we previously identified to be highly variable in both natural and experimental C. elegans populations ([36] and see Additional file 10 for details).

Data storage

Microarray data (both RNA and DNA) can be found at [23].