Genetic Variation in Parthenogenetic Collembolans Is Associated with

Dec 20, 2012 - variation in parthenogenetic populations of F. candida, and (2) this variation affects ... of females.23 Parthenogenesis is most likely...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/est

Genetic Variation in Parthenogenetic Collembolans Is Associated with Differences in Fitness and Cadmium-Induced Transcriptome Responses Benjamin Nota,*,†,‡ Maarten de Korte,† Bauke Ylstra,§ Nico M. van Straalen,† and Dick Roelofs† †

Department of Animal Ecology, Institute of Ecological Science, VU University, Amsterdam, The Netherlands Department of Clinical Chemistry, VU University Medical Center, Amsterdam, The Netherlands § Department of Pathology, VU University Medical Center, Amsterdam, The Netherlands ‡

S Supporting Information *

ABSTRACT: Ecotoxicological tests may be biased by the use of laboratory strains that usually contain very limited genetic diversity. It is therefore essential to study how genetic variation influences stress tolerance relevant for toxicity outcomes. To that end we studied sensitivity to cadmium in two distinct genotypes of the parthogenetic soil ecotoxicological model organism Folsomia candida. Clonal lines of both genotypes (TO1 and TO2) showed divergent fitness responses to cadmium exposure; TO2 reproduction was 20% less affected by cadmium. Statistical analyses revealed significant differences between the cadmium-affected transcriptomes: i) the number of genes affected by cadmium in TO2 was only minor (∼22%) compared to TO1; ii) 97 genes showed a genotype × cadmium interaction and their response to cadmium showed globally larger fold changes in TO1 when compared to TO2; iii) the interaction genes showed a concerted manner of expression in TO1, while a less coordinated pattern was observed in TO2. We conclude that (1) there is genetic variation in parthenogenetic populations of F. candida, and (2) this variation affects life-history and molecular end points relative to cadmium toxicity. This sheds new light on the sources of biological variability in test results, even when the test organisms are thought to be genetically homogeneous because of their parthenogenetic reproduction.



INTRODUCTION A major concern in ecotoxicology is that genetic variation among populations of test species can influence the outcome of ecotoxicity tests.1,2 The hypothesis that laboratory population responses are different from responses from field populations is supported by several lines of evidence from invertebrates3−5 to vertebrates.6 The aquatic ecotoxicological model Daphnia magna, for instance, shows considerable genetic variation in cadmium tolerance within and among populations. Barata et al. showed that tolerant phenotypes were not associated with costs regarding fitness,3 implying that laboratory strains grown under optimal conditions will show similar or lower tolerance to this stressor. It was suggested that this would cause an increased uncertainty in extrapolating lab results to more realistic stress scenarios in the field and that genetic variation in toxicant susceptibility should be considered more explicitly in ecological risk assessment. In contrast, survival and reproduction of different clones (or strains) of the parthenogenetic springtail species Folsomia candida Willem 1902 showed only small interclonal variation.7,8 It was concluded that the differences found were not large enough to influence the outcome of standardized soil ecotoxicity testing, such as described in ISO or OECD protocols.9,10 Thus it seems difficult to quantify the uncertainty associated with extrapolation of stress responses from the laboratory to the field. © 2012 American Chemical Society

Thus far, most of the studies considering genetic variation in ecotoxicology have focused on variation in life-history end points at the phenotypic level. In such cases, genetic variation is assumed to be present but has not been quantified at the molecular level. This is surprising, because genetic and genomic tools are proposed to strengthen ecotoxicological research11 and have subsequently been successfully implemented for certain model laboratory strains.12,13 Especially, gene expression profiling (transcriptomics) with microarrays is a powerful technique that can identify toxicity mechanisms of environmental soil pollutants.14−17 In the case of genetic adaptation of the springtail Orchesella cincta, we have shown that natural populations differ in tolerance to metals.18 Genetic variation in the metallothionein promoter region was associated with metal tolerance.19,20 Transcriptomic analyses revealed that metaltolerant springtails showed a constitutive pattern of overexpression for particular genes regardless of the level of cadmium exposure, compared to nontolerant animals. Furthermore, the overall transcriptome of tolerant animals was much less affected by metal toxicity,21 suggesting an Received: Revised: Accepted: Published: 1155

October 1, 2012 December 20, 2012 December 20, 2012 December 20, 2012 dx.doi.org/10.1021/es303983z | Environ. Sci. Technol. 2013, 47, 1155−1162

Environmental Science & Technology

Article

combined in contigs using EMBOSS 6.3.1: merger (http:// mobyle.pasteur.fr/cgi-bin/portal.py?#forms::merger), using default settings. Contigs of the different clones were aligned using ClustalW 2.0.12 (http://mobyle.pasteur.fr/cgi-bin/portal. py?#forms::clustalw-multialign), using default settings. Sequencing of other genes of interest relative to cadmium tolerance were Fcc00863, Fcc01639, Fcc02623, methallothionein-like gene (mtc), and filamenting temperature-sensitive mutant Z gene (FtsZ) and was performed as described above, and primers of all PCRs are summarized in Table S1. All contig sequences are deposited to GenBank (JQ250774-JQ250794). Toxicity Tests. Toxicity experiments were performed following standardized protocol9 and as described previously.31 In short, ten age-synchronized animals (from TO1 and TO2) were exposed to either control (Cd-) LUFA 2.2 soil (Speyer, Germany) or to cadmium containing (Cd+) LUFA 2.2 soil (spiked with 41.2 mg cadmium kg−1 soil, in the form of Cd(NO3)2·4H2O). This concentration is equivalent to the 50% effective concentration (EC50) for inhibiting reproduction in clone DK.32 Five biological replicates were used, each consisting of 10 age-synchronized animals added to a glass jar filled with 30 g of soil. After 28 days, the number of produced juveniles was determined for each jar. Concentrations of cadmium were measured in control as well as spiked soil using nitric acid extraction and atomic absorption spectrophotometry (AAS), as described previously.33 Microarray Experiment. Microarray analysis was performed to measure transcriptomic variation between the two genotypes in TO. For microarray analysis; animals of both TO1 and TO2 were exposed to “Cd-” and “Cd+” LUFA 2.2 soil (concentration cadmium see above), and total RNA was isolated according to Nota et al.32 In short, each clone/ treatment condition consisted of four replicate jars, each containing approximately 30 age-synchronized animals (23 days old). The animals were exposed for two days at the same experimental conditions as the toxicity experiment. Previous studies showed many genes to respond to cadmium exposure after a two day exposure.16 All animals per jar were extracted from soil and snap frozen in liquid nitrogen, and RNA was extracted, which included a DNase treatment. Approximately 300 ng input of total RNA was used for amplification and labeling with the Low Input Quick Amp Labeling Kit, two-color (Agilent Technologies), according to the manufacturer’s guidelines. Labeled and amplified cRNA was purified, hybridized, washed, and scanned as described in Nota et al.32 The hybridization design is shown in Figure S1. Custom Gene Expression Microarrays (Agilent Technologies) with 8 × 15 k format were used, which were earlier described in Nota et al.31 This array contains 5,069 features printed in triplicate, representing a unique gene contig from Collembase (http:// www.collembase.org/).30 More detailed description of the design of this microarray (platform) is available from Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) under accession number GPL7150. Microarray Data Analysis. Spot intensities were quantified with Feature Extraction (9.5.1.1) software (Agilent Technologies), and an overall impression of the quality of data was examined by making MA-plots. Furthermore, the “Expected LogRatio vs Observed LogRatio” of the Agilent spike-in control probes were examined and showed acceptable R2 values between 0.978 and 0.993. Preprocessing and normalization were done in the “limma” (3.6.6) package from R (2.12.0).34 This consisted of background correction (normexp, offset =

important role for transcriptional regulation in evolution of stress tolerance. In this study we investigated the effect of genetic variation upon life-history and molecular end points, within a parthenogenetic strain of F. candida: TO.22 This springtail species is an important soil ecotoxicological model and considered to be an obligate parthenogen, consisting exclusively of females.23 Parthenogenesis is most likely imposed by Wolbachia infection and is maintained by the fact that the presence of Wolbachia is required for eggs to develop.24 However, sequenced nuclear and mitochondrial markers have demonstrated genetic variation in TO, despite its assumed monoclonal genetic background. We compared the two genetically divergent TO lines against the monoclonal strain (DK) most extensively used in ecotoxicology, as this strain was genetically characterized previously.25−29 To have more insight in the different genotypes, the effects and consequences of the different genetic make-ups within the TO strain are investigated using life-history and molecular end points. Our data provide strong evidence that differences in genetic background do exist and are associated with differences in reproduction success under stressful conditions for which we also found signatures at the transcriptional level.



MATERIALS AND METHODS Animal Rearing and Genotyping. Strains of F. candida from France (TO) and Germany (DK) were established in our laboratory in 2004. Their origin is earlier described by Tully et al.22 Both have been maintained at 15 °C in 12:12 h light/dark cycles and were fed baker’s yeast, according to Nota et al.,16 following standardized methods.9 In 2007 new clonal lines were established and genotyped to ensure the purity of the clones in our experiments. This was done by selecting (5−8) individuals randomly from each strain (stock culture); each individual was kept separately in plastic rings (with lids) containing moisturized plaster of Paris mixed with charcoal, at 20 °C. Most isolated individuals deposited egg clutches and therefore started new clonal lines, which were genotyped by partially sequencing of two genes: the mitochondrial gene COII (Collembase30 (http://www.collembase.org/) accession number: Fcc01457) and gene Fcc00584. The sequence of Fcc00584 has no significant blast annotation; therefore we presume that this is not a mitochondrial but a nuclear gene. These two genes showed variation between different strains of F. candida in previous unpublished studies. DNA was isolated from 2−3 pooled animals from the newly formed clonal lines, to achieve concentrations suitable for Sanger sequencing, using the Wizard SV Genomic DNA purification kit (Promega). DNA isolation procedure was performed again just before the toxicity experiment, as a second verification of the purity of the clonal lines. Polymerase chain reaction (PCR) was performed to amplify the two genes using a mixture of MRC taq polymerase (MRC-Holland; 1 U per reaction) and Pfu DNA polymerase (Promega; 0.1 U per reaction). Primers are summarized in Table S1. PCR-products, separated by electrophoresis, were cut out of agarose gel and purified using Wizard SV Gel and PCR Clean-Up System (Promega) according to manufacturer’s protocol. 10−100 ng of PCR-product was used as template for Big Dye Cycle Sequencing (v1.1; Applied Biosystems) on an ABI PRISM 3100-Avant Genetic Analyzer. Sequence quality was manually assessed in Chromas Lite 2.0 (Technelysium Pty Ltd.), and primer sequences or low quality ends were trimmed off. Forward and reverse sequences of the each clone were 1156

dx.doi.org/10.1021/es303983z | Environ. Sci. Technol. 2013, 47, 1155−1162

Environmental Science & Technology

Article

ANOVA revealed a significant effect of the cadmium treatment and a significant genotype × treatment interaction (F = 6.12, p = 0.025, Table 1).

50), loess normalization, and quantile normalization between the arrays. Different contrasts were made, for each contrast an average (log2 transformed) ratio and a Benjamini-Hochberg adjusted p-value was generated per gene (adjusted p < 0.05 was considered significant). Limma’s RG.MA() function was used to extract the normalized intensities of the interaction genes of each hybridized sample (channel). The intensities were log2 transformed and used to generate correlation plots for each of the clones separately. These normalized (log2) intensities were also used for principal component analysis (PCA), as earlier described in Hawliczek et al.,35 using FactoMineR (1.14) in R. Annotation of the genes was performed in Blast2GO as described by Janssens et al.25 Raw and processed microarray data are available from Gene Expression Omnibus (http:// www.ncbi.nlm.nih.gov/geo/) under accession number GSE34167. Reverse Transcription and Quantitative PCR. Quantitative reverse transcription PCR (qPCR) for mtc was performed according to Nota et al.,33 using cDNA synthesized from the same RNA samples as used for microarray analysis. One reference gene (YWHAZ; Fcc02512) was used for normalization according to Nota et al.16 To verify again the purity of the TO clonal lines that were used for microarray analysis and qPCR, the nuclear marker Fcc00584 was sequenced from the cDNA samples. Two cDNA samples per clonal line were genotyped for this verification and confirmed that the clonal lines still represented the two previously identified pure TO genotypes. Statistical Analysis. For the toxicity experiment two-way analysis of variance (ANOVA) was used with clone (genotype) and treatment (cadmium) as main effects. A two-way ANOVA was used to determine the effects of cadmium upon expression of mtc, and mean normalized expression (MNE) values were first log2 transformed in order to get equal variance and normal distributed data. Levene’s test was used to assess equal variance between the treatments/genotypes and the Shapiro Wilk’s test on residuals was performed to assess normality of the data.36 All statistical analyses were performed using R (2.12.0).

Figure 1. Norms of reaction for TO1 and TO2 in control (Cd-) or cadmium (Cd+) condition. On the y-axis, the mean number of juveniles produced after 28 days (n = 5) and on the x-axis the treatment (condition). Error bars indicate standard errors.

Table 1. Results of Two-Way ANOVA for Effects of Genotype and Cadmium on Juvenile Production in Folsomia candida in a 28 Days Soil Test condition genotype cadmium treatment genotype × cadmium interaction residuals



RESULTS Influence of Genotype and Cadmium on Reproduction. Individuals of the French (TO) and the Berlin (DK) strains were isolated in order to generate pure clonal lines. Two marker genes (mitochondrial gene COII and Fcc00584) were used to genotype the two strains and showed genetic variation within TO. All individuals within TO (i.e., clonal lines) had an identical genotype (TO1) except for one (TO2). Both TO genotypes were different from the only genotype found in DK (Figures S2 and S3). The sequence of the mitochondrial gene COII was identical in TO1 and TO2, which differed from DK. However, the sequence of Fcc00584 was identical between DK and TO2 and differed from TO1 by two insertions and one point deletion. Variation in reproduction was determined relative to genotype with regard to metal toxicity by exposing them to control (Cd-) and cadmium containing (Cd+) soil for 28 days. The nominal cadmium concentration that was used is 41.2 mg kg−1 soil, which is equivalent to the EC50 on reproduction after 28 days, estimated with the DK strain.32 Concentrations of cadmium determined by AAS were 0.09 mg kg−1 in control soil and 41.18 mg kg−1 in spiked soil. Cadmium also caused a reduction in reproduction of 50% in TO1, but only of 30% in TO2, compared to the control soil (Figure 1). A two-way

Df

MS

F-value

p-value

1 1 1

570544 16820 64525

1.5976 54.1902 6.1286

0.22435 1.604 × 10−6 0.02486

16

10529

Influence of Genotype and Cadmium on Transcriptome. Microarrays were used to measure transcript levels of 5,069 genes of F. candida in TO1 and TO2, after two days of exposure to either “Cd-” or “Cd+” soil. Contrasts were made for different conditions to examine significant differential gene expression at an adjusted false discovery rate (FDR) of p < 0.05 (summarized in Table 2). There is a remarkably large difference in the number of significantly regulated genes between the two genotypes in response to cadmium exposure (1025 in TO1 vs 223 in TO2). In accordance with the smaller effects of cadmium on reproduction, the transcriptome of TO2 seems also much less affected by cadmium. Furthermore, gene expression among a number of genes was found to vary relative to genotype among control and cadmium-exposed animals. The lists of differentially expressed genes associated with the main terms cadmium and genotype can be found in Table S2−8, and the number of overlapping genes are presented in Venn diagrams in Figure S4. The expression of in total 97 genes showed a significant genotype × treatment interaction, indicating that the cadmiumaffected expression of these genes was dependent on the genotypic background. Figure 2 shows a dot plot of the log2 fold changes in response to cadmium of these interaction genes 1157

dx.doi.org/10.1021/es303983z | Environ. Sci. Technol. 2013, 47, 1155−1162

Environmental Science & Technology

Article

did not show a clear pattern of correlated gene expressions (Figure 3B). To further investigate coregulated gene expression, we performed principal component analysis (PCA) on the 97 interaction genes. Principal components (dimensions) 1, 2, and 3 explained 53.3%, 11.7%, and 8.5% of the data, respectively. The genes that are significantly correlated with principal components 1, 2, and/or 3 are summarized in Table S9. Cadmium treatment and genotype were significantly correlated with principal components 1 and 3, respectively (Figure S5, Table S9). Genetic Bases of Transcriptional Variation. We tested whether the differences of expression of the 97 interaction genes were mainly caused by transcriptional regulation, or whether sequence variation within the coding region was also involved. Therefore, we sequenced fragments of three “interaction” genes (all nuclear); a fasciclin II (Fcc00863), a monooxygenase DBH-like 1 (Fcc01639), and an endophilin A (Fcc02623). Since the microarray used in this study was designed with available sequence information from the DK clone,16,30 we also sequenced the genes from this strain. Alignments of these genes are shown in Figures S6−8. The sequences of fasciclin II and endophilin A are identical across all three clones, but sequence variation could be observed in monooxygenase DBH-like 1. The variation in this gene supports genetic variation observed for Fcc00584, confirming that TO2 and DK are identical with regard to nuclear DNA, and differed by 37 single point mutations and a 4 bp deletion from TO1. Furthermore, we partially sequenced the gene encoding the deduced metallothionein-like motif containing protein (MTC), which is (highly) upregulated in response to cadmium,37 and many other metals.33 A part of the promoter and coding region was sequenced for all three clonal lines, and the results are shown in Figure S9. In contrast to genetic variation identified in the monooxygenase and Fcc00584, mtc of the two TO clones are identical but different from DK. This suggests that not all nuclear markers are identical between DK and TO2, implying mixing of both (DK and TO1) genotypes. Finally, we sequenced Wolbachia’s FtsZ gene partially (Figure S10). Since Wolbachia is an endosymbiont and involved in F. candida’s parthenogenetic reproduction, it is interesting to examine the variation (and origin) of this bacteria in the three clonal lines. It appears that the TO clones share the same Wolbachia strain, which is different from the one previously found in DK. Expression of mtc. Due to the absence of a methallothionein-like gene (mtc) probe on our array platform the expression of this gene could not be measured in the transcriptomic study. Since methallothionein plays such an important role in cadmium resistance,20,21 we used qPCR to quantify mtc expression in the RNA samples from the microarray experiment. The mean normalized expression (MNE) values are shown in Figure S11. In response to cadmium mtc was 419.6 and 365.3 fold upregulated in TO1 and TO2, respectively. The background level of mtc (hence in “Cd-“ condition) was 2.8 fold higher in TO2 compared to TO1. Two-way ANOVA revealed a significant effect of the cadmium treatment and genotype but no significant genotype × treatment interaction effect (Table 3). A previous study with similar experimental conditions using DK showed a 602.9 fold upregulation of mtc in response to cadmium.33

Table 2. Number of Significantly Differentially Expressed Genes (FDR-Adjusted p < 0.05) in Different F. candida Strains (TO1 and TO2) Exposed to Cadmium condition cadmium: overall effect cadmium: effect in TO1 cadmium: effect in TO2 genotype: overall effect genotype: effect in Cdgenotype: effect in Cd+ interaction

contrast

total

up

down

(Cd+)-(Cd-)

1050

457

593

(Cd+)-(Cd-)

1025

438

587

(Cd+)-(Cd-)

223

129

94

(TO1)-(TO2)

991

519

472

(TO1)-(TO2)

579

269

310

(TO1)-(TO2)

668

368

300

(TO2[Cd+]-TO2[Cd])-(TO1[Cd+]-TO1[Cd-])

97

83

14

Figure 2. Dot plot of the expression (M = log2 fold change) of the 97 interaction genes of TO1 (blue) and TO2 (pink). Note, the genes were ordered from high to low using TO2’s M-values. Accession numbers of the genes are listed left, derived from Collembase.

for both genotypes. Clearly, the magnitude of fold change is in general larger in TO1 genotype as compared to TO2. To examine if the interaction genes were coregulated, we calculated the correlations of gene expressions between the 97 interaction genes for both genotypes separately and displayed them in correlation plots (Figure 3). The correlation plot for TO1 indicates a strong pattern of positively and negatively correlated gene expressions (Figure 3A), suggesting clustering of genes in coregulated networks. In contrast, the plot for TO2 1158

dx.doi.org/10.1021/es303983z | Environ. Sci. Technol. 2013, 47, 1155−1162

Environmental Science & Technology

Article

Figure 3. Correlation plots of the 97 interaction genes. Correlations (Pearson) between the expressions of two genes across all biological samples of (A) TO1 and (B) TO2; red boxes are positive correlations, blue is negative (anticorrelation), white is no correlation. Each matrix is symmetric through the diagonal, and therefore the gene order from top to down is the same from the order from left to right. Notice that the gene order in A is different from the order in B, since the order is dependent on hierarchical clustering.

(correlation) in TO1 but not in TO2. This indicates that the transcriptional response to cadmium in TO1 is in a concerted manner and highly responsive as a result of the adverse effect caused by cadmium. In contrast, TO2 is likely able to maintain homeostasis under the same cadmium toxic pressure. Globally, the processes invoked by cadmium in both TO genotypes revealed by the annotation/functions of the differentially expressed genes (Tables S3−4) are comparable. The Venn diagram in Figure S4A illustrates also that most of the cadmium affected genes in TO2 overlap with those in TO1. We find genes involved (among others) in redox reactions, oxidative stress, and antibiotics production. This is in accordance with earlier described transcriptional responses to cadmium in clone DK.16 The key differences between TO1 and TO2 are best illustrated by the differences in responsiveness to cadmium among the interaction genes shown in the dot plot (Figure 2). The fold changes are mostly much higher in TO1, which seems therefore more responsive to cadmium when compared to cadmium-exposed TO2. Moreover, the fold change of mtc regulation in response to cadmium is higher in TO1 compared to TO2, although the absolute transcript level seems to be higher for TO2 in control condition as well as in cadmium-exposed condition. The promoter and partial coding sequence of mtc did not show genetic variation among the TO clones, suggesting that the differences in mtc transcript abundance between the TO strains is caused by trans-acting factors. However, the presence of cis-regulatory variation upstream of the sequenced promoter region cannot completely be excluded and may influence the observed variation in mtc transcriptional activity. Furthermore, a constitutive high level of cadmium-protective genes (hence in “Cd-” condition), such as the metallothionein-like gene, could explain lower fold changes and was earlier postulated as a protective mechanism in the cadmium-tolerant springtail populations of O. cincta.21 Such a mechanism might cause a trade-off in the total number of TO2 juveniles produced under normal conditions, since constitutive expression of genes will cost extra energy. PCA is an efficient way of reducing complexity in transcriptomic data. The principal components (dimensions) 1 and 3 were significantly correlated with cadmium treatment

Table 3. Results of Two-Way ANOVA for Effects of the Genotype and Cadmium on mtc Expression in Folsomia candida condition genotype cadmium treatment genotype × cadmium interaction residuals

Df

MS

F-value

p-value

1 1 1

6.107 288.335 0.376

13.5638 640.4471 0.8355

0.003132 8.802 × 10−12 0.378689

12

0.450



DISCUSSION Our results show that a parthenogenetic strain TO of F. candida unexpectedly contained at least two genotypes that were associated with different levels of cadmium tolerance. The global gene expression pattern of the most cadmium sensitive genotype (TO1) showed strong signals of stress-induced transcriptional perturbation, while the tolerant genotype (TO2) retained homeostasis at the transcriptional level despite cadmium exposure The significant genotype × treatment interaction suggests that the effect of cadmium on reproduction is dependent on genotype. From the response curves (Figure 1) it is clear that the slopes have different values and cross each other. The response of TO1 shows a very steep slope, indicating a highly plastic reaction to cadmium. TO2’s response is less steep and thus less affected by cadmium. However there seems to be a trade-off in performance in control (Cd-) and treatment (Cd+) situation (i.e., TO2 has less juveniles in “Cd-” but more in “Cd +” compared to TO1). The transcriptomic response to cadmium shows even more pronounced differences between the TO clones. First, the large number of genes significantly responsive to cadmium in TO1 indicates that cadmium is more toxic to this genotype, which matches with the reaction norm at the level of reproduction (Figure 1). Second, the interaction genes are in general less responsive to cadmium (i.e., lower fold changes) in TO2 than in TO1 (Figure 2), suggesting that stress-inducibility is decreased in the more cadmium tolerant TO2 strain. Third, the interaction genes show a clear pattern of coregulation 1159

dx.doi.org/10.1021/es303983z | Environ. Sci. Technol. 2013, 47, 1155−1162

Environmental Science & Technology

Article

and genotype, respectively. Therefore the expression of the genes significantly correlated with dimension 1 is mostly involved in the cadmium response. Again here we find genes involved in redox reactions (monooxygenases, retinol dehydrogenases) but also cuticular proteins. The cuticular proteins are most likely involved in molting of (gut) epithelium, which has been shown to be an effective mechanism for excretion of cadmium from springtails.21 Genes correlated with dimension 3, and thus with genotype, are fasciclin II and apolipophorin, which are involved in cell adhesion and fat metabolism, respectively. Their exact role is unclear. When functions of the 28 annotated interacting genes are compared to the interacting genes involved in cadmium tolerance in the sister species O. cincta, striking similarities are observed. Cuticle formation, signaling via the stressactivated protein kinase (SAPK) pathway, regulation of redox state, and stress response are biological processes clearly affected in O. cincta.21 Here we identified two genes involved in cuticle formation (cuticle protein isoform b and cuticle protein 2), three genes involved in SAPK signaling (cAMP-dependent protein kinase type II, tyrosine kinase, and tyrosine-protein kinase transmembrane receptor), five redox state-related genes (thioredoxin, short chain dehydrogenase reductase, and three monooxygenases), and the stress responsive mtc gene. Although the F. candida TO strain did not evolve metal tolerance in the environment, it exerts similar regulatory changes associated with cadmium tolerance when compared to O. cincta populations that evolved cadmium tolerance at metal contaminated sites. This emphasizes the importance of such commonalities across species with regard to heavy metal detoxification. An explanation for the observed genetic variation remains to be elucidated. We showed that the sequence of the mitochondrial COII gene and Wolbachia’s FtsZ gene from the two TO clones were identical, which implies that the maternal origin is from TO. It is interesting that the Wolbachia strain is still the same in both TO clones and different from DK, since Wolbachia infection is associated with, and necessary for, parthenogenesis.24,38 Sexual reproduction may have caused generation of genetic variation, but this can only occur in the absence of Wolbachia infection. If this would have occurred in TO we expect its offspring to be reinfected with Wolbachia in the animal culture, possibly via horizontal transfer just like described in wasps.39 Males are a rare phenomenon within F. candida; we observed males very rarely, usually when cultures suffered from mites or (harmful) fungal infections, suggesting that biotic stress can trigger escape from Wolbachia and generation of males. This is also observed in Daphnia, where biotic and abiotic stress are important triggers for switching from parthenogenetic to sexual reproduction.40 In conclusion, our data show that parthenogenetic strains of F. candida can contain different genotypes. Furthermore we show that genetic variation influences the outcome of ecotoxicological tests and also of transcriptomics experiments. Previous studies with different strains of F. candida revealed small interclonal differences in response to toxicants.7,8,41 However, our results show that differences within one strain significantly alter the outcome of ecotoxicity tests. This is the first study providing a mechanistic explanation for these interclonal differences: stress-tolerant phenotypes show a constitutive overexpression of genes involved in protective pathways. Such a mechanism was already proposed with regard to populations genetically adapted to stressful environments.42

The fundamental difference, however, is that the TO2 clone did not originate from a contaminated site,22 where adaptive evolution of stress tolerance could have occurred. It implies that stress resistance can be obtained in genetically different clonal lines without the need for a stress factor as selective agent. This can in turn affect ecotoxicological end points, and it is therefore recommended that regular genetic testing of laboratory strains should take place in ecotoxicological research centers to avoid spurious outcome. Finally, the challenge to integrate data, as presented in this study, into future ecological risk assessment remains. We propose that a systems biology approach using computational models will become essential to achieve this. Such an approach was pioneered by Williams et al. using the European Flounder,43 which is a model for coastal water quality assessment. They were indeed able to predict health impacts of environmental pollution by constructing network modules linking transcriptional and metabolic responses to the physiology of compensatory adaptation in wild flounder populations. This may provide the opportunity to study environmental adaptation in a wide range of hazard scenarios in wild populations.



ASSOCIATED CONTENT

S Supporting Information *

This study contains additional figures and tables. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions

B.N. and M.d.K. performed the toxicity and microarray experiments. B.N. performed all other experiments and analyses. All authors were involved in the general setup and design of this research. D.R. supervised the project. B.N. wrote the manuscript with contributions from D.R., N.v.S., and B.Y. All authors read and approved the final manuscript. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We would like to thank Janine Mariën (VU University Amsterdam) and Thierry Janssens (BioDetection Systems Amsterdam) for bioinformatic support; Elaine van Ommen Kloeke and Bertanne Visser (VU University Amsterdam) for statistical advice; Thomas Tully (IUFM Paris) for the information on the strains; and Tjalf de Boer (VU University Amsterdam) for presenting this study at SETAC North America. Special thanks to Juriaan Buikema (Buikema Muziek Utrecht) for correct preparation of the figures.



REFERENCES

(1) Baird, D. J.; Barber, I.; Bradley, M.; Soares, A. M. V. M; Calow, P. A comparative study of genotype sensitivity to acute toxic stress using clones of Daphnia magna straus. Ecotoxicol. Environ. Saf. 1991, 21, 257−265. (2) Bickham, J. The four cornerstones of evolutionary toxicology. Ecotoxicology 2011, 20, 497−502. (3) Barata, C.; Baird, D. J.; Mitchell, S. E.; Soares, A. M. V. M. Among- and within-population variability in tolerance to cadmium stress in natural populations of Daphnia magna: Implications for ecological risk assessment. Environ. Toxicol. Chem. 2002, 21, 1058− 1064.

1160

dx.doi.org/10.1021/es303983z | Environ. Sci. Technol. 2013, 47, 1155−1162

Environmental Science & Technology

Article

(4) Agra, A.; Soares, A. V. M; Barata, C. Life-history consequences of adaptation to pollution. “Daphnia longispina clones historically exposed to copper”. Ecotoxicology 2011, 20, 552−562. (5) Timmermans, M. J. T. N.; Ellers, J.; Roelofs, D.; Van Straalen, N. M. Metallothionein mRNA expression and cadmium tolerance in metal-stressed and reference populations of the springtail Orchesella cincta. Ecotoxicology 2005, 14, 727−739. (6) Marquis, O.; Miaud, C.; Ficetola, G. F.; Bocher, A.; Mouchet, F.; Guittonneau, S.; Devaux, A. Variation in genotoxic stress tolerance among frog populations exposed to UV and pollutant gradients. Aquat. Toxicol. 2009, 95, 152−161. (7) Crommentuijn, T.; Stab, J. A.; Doornekamp, A.; Estoppey, O.; Vangestel, C. A. M. Comparative ecotoxicity of cadmium, chlorpyrifos and triphenyltin hydroxide for 4 clones of the parthenogenetic collembolan Folsomia-candida in an artificial soil. Funct. Ecol. 1995, 9, 734−742. (8) Chenon, P.; Rousset, A.; Crouau, Y. Genetic polymorphism in nine clones of a parthenogenetic collembolan used in ecotoxicological testing. Appl. Soil Ecol. 2000, 14, 103−110. (9) ISO Soil quality - Inhibition of reproduction of Collembola (Folsomia candida) by soil pollutants, 1999. (10) OECD Test No. 232: Collembolan Reproduction Test in Soil; OECD Publishing: Paris, 2009. (11) Ankley, G.; Daston, G.; Degitz, S.; Denslow, N.; Hoke, R.; Kennedy, S.; Miracle, A.; Perkins, E.; Snape, J.; Tillit, D.; Tyler, C.; Versteeg, D. Toxicogenomics in regulatory ecotoxicology. Environ. Sci. Technol. 2006, 40, 4055−4065. (12) Van Straalen, N. M.; Feder, M. E. Ecological and evolutionary functional genomicshow can it contribute to the risk assessment of chemicals? Environ. Sci. Technol. 2011, 46, 3−9. (13) Brulle, F.; Morgan, A. J.; Cocquerelle, C.; Vandenbulcke, F. Transcriptomic underpinning of toxicant-mediated physiological function alterations in three terrestrial invertebrate taxa: A review. Environ. Pollut. 2010, 158, 2793−2808. (14) Gong, P.; Guan, X.; Inouye, L.; Pirooznia, M.; Indest, K.; Athow, R.; Deng, Y.; Perkins, E. Toxicogenomic analysis provides new insights into molecular mechanisms of 2,4,6-trinitrotoluene in Eisenia fetida. Environ. Sci. Technol. 2007, 41, 8195−8202. (15) Owen, J.; Hedley, B.; Svendsen, C.; Wren, J.; Jonker, M.; Hankard, P.; Lister, L.; Sturzenbaum, S.; Morgan, A.; Spurgeon, D.; Blaxter, M.; Kille, P. Transcriptome profiling of developmental and xenobiotic responses in a keystone animal, the oligochaete annelid Lumbricus rubellus. BMC Genomics 2008, 9, 266. (16) Nota, B.; Timmermans, M. J. T. N.; Franken, O.; MontagneWajer, K.; Mariën, J.; De Boer, M. E.; De Boer, T. E.; Ylstra, B.; Van Straalen, N. M.; Roelofs, D. Gene expression analysis of Collembola in cadmium containing soil. Environ. Sci. Technol. 2008, 42, 8152−8157. (17) Amorim, M. J. B.; Novais, S. C.; Van Der Ven, K.; Vandenbrouck, T.; Soares, A. M. V. M.; De Coen, W. Development of a microarray for Enchytraeus albidus (Oligochaeta): preliminary tool with diverse applications. Environ. Toxicol. Chem. 2011, 30, 1395− 1402. (18) Janssens, T. K. S.; Roelofs, D.; Van Straalen, N. M. Molecular mechanisms of heavy metal tolerance and evolution in invertebrates. Insect Sci. 2009, 16, 3−18. (19) Janssens, T. K. S.; Lopez, R. D. R.; Marien, J.; Timmermans, M. J. T. N.; Montagne-Wajer, K.; Van Straalen, N. M.; Roelofs, D. Comparative population analysis of metallothionein promoter alleles suggests stress-induced microevolution in the field. Environ. Sci. Technol. 2008, 42, 3873−3878. (20) Janssens, T. K. S.; Marien, J.; Cenijn, P.; Legler, J.; Van Straalen, N. M.; Roelofs, D. Recombinational micro-evolution of functionally different metallothionein promoter alleles from Orchesella cincta. BMC Evol. Biol. 2007, 7, 88. (21) Roelofs, D.; Janssens, T. K. S.; Timmermans, M. J. T. N.; Nota, B.; Marien, J.; Bochdanovits, Z.; Ylstra, B.; Van Straalen, N. M. Adaptive differences in gene expression associated with heavy metal tolerance in the soil arthropod Orchesella cincta. Mol. Ecol. 2009, 18, 3227−3239.

(22) Tully, T.; D’Haese, C. A.; Richard, M.; Ferriere, R. Two major evolutionary lineages revealed by molecular phylogeny in the parthenogenetic Collembola species Folsomia candida. Pedobiologia 2006, 50, 95−104. (23) Fountain, M. T.; Hopkin, S. P. Folsomia candida (Collembola): A “standard” soil arthropod. Annu. Rev. Entomol. 2005, 50, 201−222. (24) Timmermans, M. J. T. N.; Ellers, J. Wolbachia endosymbiont is essential for egg hatching in a parthenogenetic arthropod. Evol. Ecol. 2009, 23, 931−942. (25) Janssens, T. K. S.; Staaden, S.; Scheu, S.; Marien, J.; Ylstra, B.; Roelofs, D. Transcriptional responses of Folsomia candida upon exposure to Aspergillus nidulans secondary metabolites in single and mixed diets. Pedobiologia 2010, 54, 45−52. (26) De Boer, T. E.; Birlutiu, A.; Bochdanovits, Z.; Timmermans, M. J. T. N.; Dijkstra, T. M. H.; Van Straalen, N. M.; Ylstra, B.; Roelofs, D. Transcriptional plasticity of a soil arthropod across different ecological conditions. Mol. Ecol. 2011, 20, 1144−1154. (27) Timmermans, M. J. T. N.; Roelofs, D.; Nota, B.; Ylstra, B.; Holmstrup, M. Sugar sweet springtails: on the transcriptional response of Folsomia candida (Collembola) to desiccation stress. Insect Mol. Biol. 2009, 18, 737−746. (28) Van Ommen Kloeke, A.; Van Gestel, C.; Styrishave, B.; Hansen, M.; Ellers, J.; Roelofs, D. Molecular and life-history effects of a natural toxin on herbivorous and non-target soil arthropods. Ecotoxicology 2012, 21, 1084−1093. (29) Nota, B.; Van Straalen, N. M.; Ylstra, B.; Roelofs, D. Gene expression microarray analysis of heat stress in the soil invertebrate Folsomia candida. Insect Mol. Biol. 2010, 19, 315−322. (30) Timmermans, M. J.; De Boer, M. E.; Nota, B.; De Boer, T. E.; Marien, J.; Klein-Lankhorst, R. M.; Van Straalen, N. M.; Roelofs, D. Collembase: a repository for springtail genomics and soil quality assessment. BMC Genomics 2007, 8, 341. (31) Nota, B.; Bosse, M.; Ylstra, B.; Van Straalen, N.; Roelofs, D. Transcriptomics reveals extensive inducible biotransformation in the soil-dwelling invertebrate Folsomia candida exposed to phenanthrene. BMC Genomics 2009, 10, 236. (32) Nota, B.; Verweij, R. A.; Molenaar, D.; Ylstra, B.; Van Straalen, N. M.; Roelofs, D. Gene expression analysis reveals a gene set discriminatory to different metals in soil. Toxicol. Sci. 2010, 115, 34− 40. (33) Nota, B.; Vooijs, R.; Van Straalen, N. M.; Roelofs, D. Expression of mtc in Folsomia candida indicative of metal pollution in soil. Environ. Pollut. 2011, 159, 1343−1347. (34) Smyth, G. K. Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 2004, 3. (35) Hawliczek, A.; Nota, B.; Cenijn, P.; Kamstra, J.; Pieterse, B.; Winter, R.; Winkens, K.; Hollert, H.; Segner, H.; Legler, J. Developmental toxicity and endocrine disrupting potency of 4azapyrene, benzo[b]fluorene and retene in the zebrafish Danio rerio. Reprod. Toxicol. 2012, 33, 213−223. (36) Shapiro, S. S.; Wilk, M. B. An analysis of variance test for normality (complete samples). Biometrika 1965, 52, 591−611. (37) Nakamori, T.; Fujimori, A.; Kinoshita, K.; Ban-nai, T.; Kubota, Y.; Yoshida, S. mRNA expression of a cadmium-responsive gene is a sensitive biomarker of cadmium exposure in the soil collembolan Folsomia candida. Environ. Pollut. 2010, 158, 1689−1695. (38) Frati, F.; Negri, I.; Fanciulli, P. P.; Pellecchia, M.; De Paola, V.; Scali, V.; Dallai, R. High levels of genetic differentiation between Wolbachia-infected and non-infected populations of Folsomia candida (Collembola, Isotomidae). Pedobiologia 2004, 48, 461−468. (39) Huigens, M. E.; Luck, R. F.; Klaassen, R. H. G.; Maas, M. F. P. M.; Timmermans, M. J. T. N.; Stouthamer, R. Infectious parthenogenesis. Nature 2000, 405, 178−179. (40) Eads, B. D.; Andrews, J.; Colbourne, J. K. Ecological genomics in Daphnia: stress responses and environmental sex determination. Heredity 2008, 100, 184−190. 1161

dx.doi.org/10.1021/es303983z | Environ. Sci. Technol. 2013, 47, 1155−1162

Environmental Science & Technology

Article

(41) Diogo, J.; Natal-da-Luz, T.; Sousa, J.; Vogt, C.; Nowak, C. Tolerance of genetically characterized Folsomia candida strains to phenmedipham exposure. J. Soils Sediments 2007, 7, 388−392. (42) Roelofs, D.; Aarts, M.; Schat, H.; Van Straalen, N. Functional ecological genomics to demonstrate general and specific responses to abiotic stress. Funct. Ecol. 2008, 22, 8−18. (43) Williams, T. D.; Turan, N.; Diab, A. M.; Wu, H.; Mackenzie, C.; Bartie, K. L.; Hrydziuszko, O.; Lyons, B. P.; Stentiford, G. D.; Herbert, J. M.; Abraham, J. K.; Katsiadaki, I.; Leaver, M. J.; Taggart, J. B.; George, S. G.; Viant, M. R.; Chipman, K. J.; Falciani, F. Towards a system level understanding of non-model organisms sampled from the environment: A network biology approach. PLoS Comput. Biol. 2011, 7, e1002126.

1162

dx.doi.org/10.1021/es303983z | Environ. Sci. Technol. 2013, 47, 1155−1162