Transcriptional plasticity of a soil arthropod across different ecological conditions

Molecular Ecology (2011) 20, 1144–1154

doi: 10.1111/j.1365-294X.2010.04985.x

Transcriptional plasticity of a soil arthropod across different ecological conditions TJALF E. DE BOER,* ADRIANA BIRLUTIU,† ZOLTAN BOCHDANOVITS,‡ MARTIJN J. T. N. TIMMERMANS,§ TJEERD M. H. DIJKSTRA,‡ NICO M. VAN STRAALEN,* BAUKE YLSTRA– and D I C K R O E L O F S * *Department of Ecological Science, VU University, de Boelelaan 1085, 1081HV, Amsterdam, the Netherlands, †Radboud University Nijmegen, Intelligent Systems, Toernooiveld 1, 6525ED, Nijmegen, the Netherlands, ‡ Department of Clinical Genetics, Section Medical Genomics, VU University Medical Center, van der Boechorststraat 7, 1081 BT, Amsterdam, the Netherlands, §Department of Ecology and Evolution, Faculty of Natural Sciences, Imperial College London, Exhibition Road, London SW7 2AZ, UK, –Department of Pathology, VU University Medical Centre, de Boelelaan 1117, 1081HV, Amsterdam, the Netherlands

Abstract Ecological functional genomics, dealing with the responses of organisms to their natural environment is confronted with a complex pattern of variation and a large number of confounding environmental factors. For gene expression studies to provide meaningful information on conditions deviating from normal, a baseline or normal operating range (NOR) response needs to be established which indicates how an organism’s transcriptome reacts to naturally varying ecological factors. Here we determine the transcriptional plasticity of a soil arthropod, Folsomia candida, exposed to various natural environments, as part of a first attempt in establishing such a NOR. Animals were exposed to 26 different field soils after which gene expression levels were measured. The main factor found to regulate gene expression was soil-type (sand or clay). Cell homeostasis and DNA replication were affected in collembolans exposed to sandy soil, indicating general stress. Multivariate analysis identified soil fertility as the main factor influencing gene expression. Regarding land-use, only forest soils showed an expression pattern deviating from the others. No significant effect of land-use, agricultural practice or soil type on fitness was observed, but arsenic concentration was negatively correlated with reproductive output. In conclusion, transcriptional responses remained within a limited range across the different land-uses but were significantly affected by soil-type. This may be caused by the contrasting soil physicochemical properties to which F. candida strongly responds. The broad range of conditions over which this soil-living detritivore is able to survive and reproduce, indicates a strategy of high plasticity, which comes with extensive gene expression regulation. Keywords: ecophysiology, genomics, natural operating range, transcriptional plasticity Received 17 May 2010; revision received 31 October 2010; accepted 3 November 2010

Introduction Genomic and gene expression measurements are becoming commonly used techniques in ecological studies (Gibson 2002; Kammenga et al. 2007). By measuring gene expression, we can determine the physiological state of animals when exposed to different Correspondence: Tjalf de Boer, Fax: +31 (0)205987123; E-mail: [email protected]

ecological conditions and determine whether some of these conditions cause adverse effects. In soil ecology one or more species of animals or plants are exposed to disturbed soils in order to determine the impact of local environmental change on soil dwelling species (Roelofs et al. 2008). To arrive at meaningful results for species functioning under stressed conditions, a reference or baseline level of functioning is needed. We dub this baseline response the normal operating range (NOR).  2011 Blackwell Publishing Ltd

S O I L A R T H R O P O D T R A N S C R I P T I O N A L P L A S T I C I T Y 1145 Similar concepts have been used in ecological studies for a long time. Odum et al. (1979) introduced the concept in their ecosystem perturbation theory already in the seventies, where they stated that perturbation is any deviation or displacement from the nominal state. The nominal state was defined not as single or fixed, but as a range of normal functioning, including expected variance, hence NOR. We propose to apply the concept of NOR to gene expression measurements. This NOR can be applied to differentiate organismal responses with negative effects on fitness, from responses to the natural environment. In complicated ecological studies, where it can be difficult to standardize environmental factors and where a large number of endpoints are measured, such a NOR will be essential to study the influence of stress factors on gene expression. Establishing a NOR however, is a daunting task. In this article, we take a first step towards the operationalization of this NOR concept by studying the plasticity of a test species with a single genotype exposed to a multitude of natural conditions. The variation of gene expression under these conditions reflects physiologically relevant plasticity and may provide an insight into how animals maintain homeostasis in the field. Transcriptome analysis has become a standard method for assessing the physiological state of an organism. Microarray technology is used to measure the expression of a large number of genes at once (Schena et al. 1995). In most microarray analyses a treatment sample (e.g. mRNA from animals subjected to some stress factor) is compared to a control sample (e.g. mRNA from a reference group). Less attention is paid to variation within control situations or reference conditions, that is, the range of gene expression variation shown by animals living under conditions that can be considered to fall within their ecological niche. Such studies will help to define baseline levels of gene expression for an organism. Pritchard et al. (2001), investigated gene expression differences among six normal male C57BL6 mice. Significant expression variation was found in 0.8–3.3% of the measured genes depending on the tissue investigated. Among the differentially expressed genes, immune-system–related genes, stressinduced genes and hormonally regulated genes were over represented. Cavalieri et al. (2000) characterized gene expression between phenotypic variants of progeny from a single parental strain and detected 6% of genes differentially expressed between these variants. Both studies emphasize the importance of determining baseline gene expression variation in order to avoid misinterpretation of microarray data. Ecological studies often focus on non-genomic model organisms, can be complex and involve a number of  2011 Blackwell Publishing Ltd

naturally varying confounding factors. Here we focus on an upcoming genomic model species with high ecological relevance for soil quality assessment, the collembolan Folsomia candida. This collembolan is an important test animal in soil ecotoxicology and is part of an ISOrecognized toxicity test (ISO 1999). F. candida is easy to culture, has a short reproduction time and has low genetic variation among individuals in the same culture due to its parthenogenetic mode of reproduction (Fountain & Hopkin 2005). At the same time, physiological differences between cultures from different origins are readily observed. These attributes make it a suitable animal for laboratory experiments on population genetics and evolution (Smit & Van Gestel 1998; Noe¨l et al. 2006). Part of the F. candida transcriptome has been sequenced and used for gene transcriptional studies (Timmermans et al. 2007; Nota et al. 2008). Timmermans et al. (2009) used an oligonucleotide microarray platform for a physiological study on the molecular mechanism of drought tolerance in this springtail. They suggested carbohydrate transport, sugar catabolism and cuticle maintenance to be important biological processes involved in combating desiccation stress. Interestingly, Bayley and Holmstrup (1999) showed that F. candida becomes hyperosmotic during desiccation by accumulating glucose and myo-inositol, so that water vapour can be extracted from the environment under desiccating conditions. Thus, the transcriptomic data supported previous physiological observations. Another study by Nota et al. (2009) investigated the biotransformation pathway of xenobiotic substances in F. candida. Indeed, genes involved in phases I, II and III of the biotransformation pathway were significantly affected by the xenobiotic compound phenanthrene. The ISO soil toxicity test measures F. candida survival and reproduction after an exposure of 28 days. This test, however, tends to show a lot of variation between replicates (Crouau & Cazes 2003) and it does not provide any information on the mode of action that a chemical or pollution stress exerts on the animal. Also, measurement of survival and reproduction does not reveal comparative information about modes of action between different test species. Measuring gene expression as an endpoint in a soil ecotoxicological test can possibly solve these issues by providing a sensitive and mode of action specific way of measuring toxicity in soil (Van Straalen & Roelofs 2008). Furthermore, gene expression analysis can be combined with survival ⁄ reproduction measurements to link molecular information to ecological endpoints. The aim of this article is to explore how gene expression in this collembolan varies under different ecological conditions that are part of its natural environment and what are the main factors determining such variation.

1146 T . E . D E B O E R E T A L . Generally, F. candida is found over a wide range of soil depth and inhabits agricultural ecosystems, forests and edges of streams (Fountain and Hopkin, 2005). Previous studies have reported on different population densities (Kaneda & Kaneko 2008) and growth rates (Kaneda & Kaneko 2002) of F. candida exposed to a variety of natural soils. To obtain more mechanistic insight into the ability of F. candida to survive and reproduce under such a broad range of conditions, we measured gene expression in animals exposed to a wide variety of natural soils. We want to contribute to the development of a baseline indicative of the transcriptional NOR in this collembolan.

Materials and methods Soil sampling Twenty-six soils from different sites in the Netherlands were sampled between March and June 2007. Soil samples were taken from fields with a known history of land-use (dairy farming, agriculture, natural grassland and forest). Information on specific agricultural practice (conventional vs. organic) for the agricultural and dairy farming land uses was also available. The soils were chosen to represent different soil types (clay vs. sand). The sites jointly represented several replicated combinations of land-use, agricultural practice and soil type, but not all combinations could be made (Table 1). F. candida, the test animal used in this study, is normally found in these kinds of soils, although we did not establish its presence in every field. Fields were sampled at five spots in a square of 20 by 20 m (on each corner and in the middle). At each spot, five subsamples were taken within a radius of 2 m (see Table 1 for soil details and map coordinates). All soils were sampled with consent of their respective owners or caretakers. Samples per soil plot were mixed, sieved over a 4 mm grid to exclude gravel and plant material and stored at 5 C. The pretreatment procedure resulted in removal of larger infauna, although the soils were not sterilized and retained their natural physicochemical properties. Chemical analysis on 250 g subsamples was performed by Blgg, Oosterbeek, the Netherlands, who measured eight metals (total Cd, Cr, Cu, Hg, Pb, Ni and As); pH; total nitrogen; total carbon; phosphate (Ptotal, Pw, P-AL and P-PAE); clay content and particle size distribution according to certified methods. We also measured the water holding capacity (WHC) for every soil. In addition to the Dutch field soils, the standard LUFA 2.2 reference soil was used (Landwirtschaftliche Untersuchungs und Forschungsanstalt, Speyer, Germany). This control soil is a sandy soil from the western part of Germany and often

used in bioassays conducted under internationally harmonized protocols.

Folsomia candida culture and exposure Folsomia candida (VU Berlin strain) was maintained in PVC containers with a plaster of Paris base containing 10% charcoal. The animals were fed baker’s yeast (Dr Oetker, Amersfoort, The Netherlands) ad libitum. Following ISO (1999). Age-synchronized cultures were obtained by transferring adult collembolans to fresh culture containers where they were allowed to lay eggs for 2 days. After 2 days the animals were removed and their hatchlings were used for the experiments. All collembolans (stocks and exposed) were kept at 20 C in a climate-controlled room with a 12-h dark ⁄ light cycle at 75% relative air humidity. Exposures were performed in 100 mL glass jars. The field soils were moistened to 50% of their WHC 24 h before exposure and left to equilibrate. For the gene expression analysis four replicate jars per soil were used. Each replicate contained 25 g wet soil and 30, 23day-old collembolans randomly selected from the synchronized stock. After 2 days of exposure, the animals were removed from the soil by floatation on water and snap frozen in liquid nitrogen for RNA extraction. The 30 animals from each jar were pooled and considered one biological replicate. For practical reasons the exposures for sandy soils and clay soils were performed consecutively. For the 28-day reproduction test, 10 day old, synchronized collembolans were used. Ten animals were exposed to 25 g of wet soil (at 50% of WHC) in a 100 mL glass container for 28 days per replicate. Containers were opened twice a week for aeration, fed once a week and moisture levels were adjusted twice during the exposure. After 28 days the containers were filled with 100 mL water and emptied in a glass beaker and after gentle stirring digital photographs were taken of the floating collembolans. The number of juveniles was determined with CellD software (Olympus, Hamburg, Germany).

RNA extraction and microarray hybridization Total RNA was extracted with the SV total RNA kit from Promega according to manufacturer’s instructions, which included a DNAse treatment. RNA integrity and concentrations were measured on a Bioanalyzer (Bioanalyzer 2100; Agilent Technologies, Santa Clara, USA) and Nanodrop spectrophotometer (Nanodrop ND-1000; Fisher Scientific, Waltham, USA). 500-ng total RNA per sample was used as input for amplification and labelling with the Low-Input Fluorescent Linear Amplification  2011 Blackwell Publishing Ltd

S O I L A R T H R O P O D T R A N S C R I P T I O N A L P L A S T I C I T Y 1147 Table 1 Specific data for each soil concerning soil-type, land-use and practice as well as GPS coordinates Soil

Soil type

Land-use

Practice

North coordinate

East coordinate

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

Sand Sand Sand Sand Sand Sand Sand Sand Sand Sand Sand Sand Clay (sea) Clay (sea) Clay (sea) Clay (sea) Clay (sea) Clay (river) Clay (river) Clay (river) Sand Sand Sand Sand Sand Clay (river)

Dairy Dairy Dairy Forest Forest Forest Agriculture Agriculture Agriculture Agriculture Agriculture Agriculture Agriculture Agriculture Agriculture Agriculture Agriculture Dairy Dairy Dairy Natural grassland Natural grassland Natural grassland Natural grassland Agriculture Dairy

Conventional Conventional Conventional N⁄A N⁄A N⁄A Conventional Conventional Conventional Organic Organic Organic Organic Organic Conventional Conventional Conventional Conventional Conventional Conventional N⁄A N⁄A N⁄A N⁄A Organic Conventional

5214¢84900 5214¢36800 5127¢99200 5208¢17300 5125¢62700 5304¢92000 5306¢33700 5243¢51500 5305¢30900 5213¢22000 5251¢71100 5201¢24100 5134¢15000 5312¢92200 5133¢00800 5252¢07600 5312¢66700 5153¢19600 5222¢86200 5129¢26800 5132¢70600 5203¢51500 N⁄A 5200¢52400 5142¢91700 5151¢26200

00616¢13000 00641¢41300 00554¢35000 00511¢18500 00547¢07400 00628¢03100 00622¢53100 00637¢16000 00649¢94500 00540¢16000 00643¢20300 00612¢04800 00335¢06600 00527¢59300 00328¢05900 00502¢06700 00531¢02100 00617¢34000 00604¢83300 00518¢11500 00518¢37500 00533¢89900 N⁄A 00535¢85300 00546¢72500 00555¢18100

N ⁄ A, not applicable.

Kit (Agilent Technologies), according to the manufacturer’s guidelines. In this protocol total RNA is used as input for reverse transcription into cDNA which in turn is used as template for labelled cRNA transcription. This results in greater amplification of the labelled material. One modification to the standard protocol was applied: the cRNA transcription reactions were done in half volume. Each 8*15K custom Agilent microarray was competitively hybridized with two colour cyanine 3 (Cy3) and cyanine 5 (Cy5) labelled cRNA samples overnight, washed and scanned, all according to manufacturer’s instructions. The custom microarray contains 5069 unique F. candida gene fragments in triplicate (Nota et al. 2009), and is based on the F. candida EST sequencing database Collembase (Timmermans et al. 2007) (http:// www.collembase.org).

Microarray experimental design and analysis An interwoven loop design (Altman & Hua 2006) was used in the microarray experiment as follows. Four replicates for each soil were labelled twice with Cy3 and twice with Cy5, obtaining a dye swap at the biological replication level. The loop design was created in such a way that replicates from one soil were never competi 2011 Blackwell Publishing Ltd

tively hybridized against more than one replicate from another soil on a single microarray. Separate loop designs were used for either clay or sandy soil exposures. In each loop a control was incorporated which consisted of collembolans exposed to the standard soil LUFA 2.2, in order to experimentally compare both exposure series. A detailed description of the design is shown in Fig. S1 (Supporting information). Microarray fluorescent intensities were measured with Feature Extraction software (version 9.5.1, Agilent technologies). The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus (Edgar et al. 2002) and are accessible through GEO Series accession number GSE21213. Pre-processing and normalization (global loess) were performed using the Limma package in R (Smyth 2004). Limma was also used for background correction (Edwards 2003) which included a minimum intensity offset of 30 to avoid zero or negative intensities. Technical replicates on the array (three per gene) were averaged before microarray ANOVA (MAANOVA) analysis for which we used the MAANOVA package in R (Kerr et al. 2000). The microarray ANOVA model was described by the following formula: yijkl = l + Si + Lj + (S x L)ij + Dk + Al + eijkl, where yijkl is the signal from the sample derived from the ith soil-type in combination with the

1148 T . E . D E B O E R E T A L . jth land-use, labelled with the kth dye and hybridized to the randomly assigned lth array. The parameter l is the overall mean, Si is the soil-type effect, Li is the land-use effect, (S x L)ij is the soil-type land-use interaction, Dk the dye effect (Cy3 or Cy5), Al is the random effect of the array to which the sample has been assigned to and eijkl is the stochastic error. In the soiltype MAANOVA results a fold regulation cut-off was determined by comparing the LUFA2.2 controls from both sand and clay exposures and this cut-off was used to remove low fold-change genes from the results. The clipped list was subjected to Gene Ontology (GO) analysis using the TopGO package in R (Alexa et al. 2006) and clustered with unsupervized hierarchical clustering (Pearson un-centred) using Tigr Mev 4.5.1 (Saeed et al. 2006). Significance cut-off of the GO terms was set at 0.05 after weighing the results from the Fisher’s exact test. Weighing was done to prevent overestimation of significance due to the hierarchically structured tree of GO terms. GO terms which, in the F. candida database were represented by only one gene were removed from further analysis. The significant gene list generated from the land-use MAANOVA calculation was subjected to K-means clustering (four groups, 100 iterations) using the MAANOVA package in R. To determine overall gene expression stability, the coefficient of variation for each gene was calculated within the soil-type ⁄ land-use combinations. Also, to determine gene expression stability between soil-type ⁄ land-use combinations a series of pair-wise F-tests were performed on variance data per soil-type ⁄ land-use subset. Calculated P-values were adjusted according to Benjamini and Hochberg’s step-up procedure (Benjamini & Hochberg 1995). A canonical correlation analysis (CCA) was performed to link gene expression to soil analytical data. CCA is a method for determining the relationship between two sets of variables and seeks linear combinations between the variables in both datasets that are maximally correlated with each other. In detail, we first factored out the sand vs. clay effect from the gene expression measurements by subtracting the mean expression level for sand or clay from each gene. Reason for doing so is that we are interested in a correlation of soil analytical data with gene expression other than the huge sand vs. clay effect. Secondly, we used a variant of penalized CCA (Witten et al. 2009) as the dimensionality of the gene expression data was larger than the number of data points and standard CCA cannot work in this case. For the analysis we used R package PMA. We used a permutation test (CCA.permute) to find the optimal L1 penalty: for penalty-x = penalty-z = 0.7 we obtained a z-statistic of 3.

The microarray platform used in this study has already been validated several times by means of QPCR in previous studies (Nota et al. 2009; Timmermans et al. 2009), resulting in a high correlation between QPCR-derived expression and microarray-derived expression (Spearman’s Rho 0.85–0.94).

Results Variation of LUFA controls in time Exposures to the two soil-types (sandy and clay soils) were performed separately in time and the samples from each exposure were contained in an interwoven loop that was used for the microarray design (Fig. S1, Supporting information). The microarray experimental design included two interwoven loops (sandy soils and clay soils). To investigate possible systematic differences between the two loops and thus the two exposures, a LUFA2.2 control was added to each exposure which was similar (moist content, pH, four replicates per exposure) in each exposure loop. These two controls were used to estimate systematic differences in gene expression potentially caused by time and different animal batches. The microarray intensities for both LUFA2.2 controls were normalized, log transformed and averaged per exposure. The two exposure controls were compared directly to each other with linear regression. This analysis revealed a slope of 0.997 and intercept of 0.077, which hardly deviates from 1 and 0 for a perfect regression where the two sets are identical. The maximum log2 fold change between the two controls was 0.42 (see Fig. 1). This shows that gene expression signatures of F. candida across standard soils are highly reproducible. This allowed for a direct comparison between clay and sand datasets in further analyses, leaving out the expression patterns from LUFA soils.

Gene expression analysis To measure gene expression differences in F. candida exposed to the different soils we determined two main levels of variation: soil-type (clay or sandy soil) and land-use (i.e. forest, organic agriculture, etc.), see Table 1 for the details per soil. A microarray ANOVA model was used to determine differences in gene expression within soil-type and land-use and to investigate if there was an interaction between these two main factors. In the main factor soil-type we found, after FDR adjustment (P < 0.05, adaptive method), 2819 out of the 5069 genes to be significantly differentially expressed between the two soil types. To exclude the possibility of an exposure effect, fold change cut-offs of 0.42 and  2011 Blackwell Publishing Ltd

S O I L A R T H R O P O D T R A N S C R I P T I O N A L P L A S T I C I T Y 1149

−2

−1

M 0

1

2

MA plot LUFA-clay vs. LUFA-sand

6

8

10

12

14

16

A

Fig. 1 An MA plot between the averages for both the clay soil exposure controls and the sandy soil exposure controls. The maximum fold-change between clay and sandy soil exposure was 0.42 (log2) which was used as a fold-change cut-off to separate transcripts regulated by soil-type from transcripts regulated by exposure effects.

)0.42 (log2) were used, since this was the maximum fold change observed between the LUFA2.2 controls of the sand and clay exposures. The number of transcripts left after the fold change cut-off was 936. These 936 genes were divided into two groups according to their response type (positive or negative) (Table S1, Supporting information). The group with a positive fold change (genes upregulated in clay soils) consisted of 449 genes and the group with a negative fold change (genes upregulated in sandy soils) contained 487 genes. Figure 2 depicts a clustering of the 936 remaining genes clearly showing contrasting profiles for clay and sandy soils. Both sets of genes were subjected to a gene enrichment (GO) analysis (Alexa et al. 2006). In the GO analysis we focused on the terms ‘biological process’ and ‘molecular function’ since these seemed to be the most indicative of any effects on a biological level. In the clay upregulated group of genes, we found for biological process GO terms that were mainly involved in protein maintenance, cytoskeleton regulation and signal transduction. Analysis on Molecular Function in the clay upregulated group revealed terms which were involved in protein stabilization, protein homeostasis and fatty acid metabolism. In the sand upregulated group both biological process and molecular function revealed a large group of GO terms involved in DNA replication ⁄ repair and the cell cycle. A smaller group which was also found in both ‘biological process’ and ‘molecular function’ included GO terms for amino acid metabolism.  2011 Blackwell Publishing Ltd

Fig. 2 Hierarchical clustering of the expression results for the 936 differentially regulated genes in soil-type. Two main clusters between exposure to sandy soils (left) and exposure to clay soils can be seen. Sand is further divided into two clusters but this division remains unexplained.

The other main factor considered in the analysis was land-use. Here we compared animals exposed to soil under organic agriculture, conventional agriculture, dairy farming, forest and natural grassland. This analysis revealed 12 genes which were significantly differentially expressed between the different land uses. K-means clustering showed specific expression patterns where forest soils caused the main effect. Four genes were found upregulated in collembolans exposed to forest soils while two genes were downregulated. According to BLAST analysis, two of the four upregulated genes were ABC transporters while the two downregulated genes are glucuronosyl transferases. Both ABC transporters and glucuronosyl transferases are part of the phases I

1150 T . E . D E B O E R E T A L . and III detoxification pathway involved in the removal of xenobiotic substances from the body. This could indicate exposure to potentially harmful organic compounds, such as polyphenols and humic acids, which are more abundant in forest soils than in agricultural soils (Hattenschwiler et al. 2005). A MAANOVA calculation was also performed where we investigated the interaction between soil-type and landuse. The analysis yielded no significant genes after FDR adjustment, indicating that the differences in gene expression between soil types are independent of landuse. It appears that the main dichotomy in the data is between clay and sand, while land-use has a smaller effect on the transcriptome, not interacting with soiltype.

Variation in gene expression

0.08

Genome-wide increased variation in gene expression may indicate stress (Oleksiak et al. 2002). We calculated the coefficient of variation (CoV) per gene for each soiltype ⁄ land-use combination in order to compare expression stability between the different soil-type ⁄ land-use combinations (Fig. 3). The CoV ranged from 0.5% to 9% with 90% of the genes showing a CoV lower than 3%. Compared to other studies which showed CoV

ranges up to 15% this is quite low (Oleksiak et al. 2002). Apparently, F. candida gene expression, as measured by the protocol we employed, is relatively stable within a single soil type, even when the soils are from different natural environments. The 30 genes with the highest CoV included two isopenicillin-N-synthase (IPNS) genes and an ACV synthase gene. In fungi and bacteria these genes are involved in antibiotic synthesis and are currently under investigation as possible novel antibiotic synthesizing genes present in soil arthropods (Nota et al. 2008). Also, two Niemann-pick type C genes were identified. In insects these genes are essential in cholesterol synthesis which is used for the production of the moulting hormone ecdyson (Huang et al. 2005). To investigate if some land uses induced greater variance in gene expression than others, a series of pair-wise F-tests was performed on soil-type ⁄ land-use combination variance data. Before the F-tests, the data was tested on normal distribution and genes that not showed a normal distribution were removed from the dataset. The F-distribution per pair-wise test was calculated for the remaining 3420 transcripts and calculated P-values were adjusted for multiple testing according to Benjamini and Hochberg’s step-up procedure. In general, low variation in gene expression was observed across the different land uses. Three pair-wise tests (forest soil vs. natural grassland on sand, conventional agriculture vs. dairy farming on sand and dairy farming on sand vs. dairy farming on clay) showed significant differences in gene expression variation (see Table 2).

Reproduction

Coefficient of variation 0.04

0.06

Reproduction of F. candida was measured in the same set of soils after a 28-day exposure. Variation in the numbers of juveniles between different soils was high. On the average, clay soils resulted in fewer juveniles Table 2 The number of significant genes in the different pair-wise tests performed to investigate the variance in gene expression variation between different soil-type ⁄ land-use combinations.

0.00

0.02

Pair-wise test

Measured transcripts

Fig. 3 Average variance per gene calculated from different soil-type ⁄ land-use combinations ordered from low variance to high variance. The majority of the genes (94%) show an average variance lower then 0.1 (10%).

Sig. genes

Sand: organic agriculture vs. conventional agriculture Sand: conventional agriculture vs. dairy farming Clay: organic agriculture vs. conventional agriculture Sand: forest vs. natural grassland Clay vs. sand organic agriculture Clay vs. sand conventional agriculture Clay vs. sand dairy farming

0 204 0 272 28 34 272

 2011 Blackwell Publishing Ltd

S O I L A R T H R O P O D T R A N S C R I P T I O N A L P L A S T I C I T Y 1151 than sandy soils (406 for clay and 558 for sand) but this was not statistically significant considering the large variation within each soil. Variation between the soils ranged from 173 juveniles on average for soil 15 (organic agriculture on clay) to 809 juveniles on average for soil 8 (conventional agriculture on sand) while the control soil LUFA2.2 contained 634 juveniles on average. We observed no statistically significant pattern in the reproduction data concerning land-use. The variation seen in the experiments (200–800 juveniles after 28 days) falls in the range which is normally observed in experiments of this type (Smit & Van Gestel 1998). Of all the environmental factors that were measured only the soil arsenic concentration was significantly correlated with the number of juveniles per soil (Pearson correlation, )0.495, P = 0.01) (Fig. 4). According to Crouau and Moı¨a (2006) the concentration of As causing 50% reduction of reproduction over 28 days (EC50) is 21.7 mg ⁄ kg soil. In one soil (soil 1, dairy farming on sand) the arsenic concentration exceeded this EC50 (26 mg ⁄ kg); this soil also had a lower reproduction (61% of the average reproduction for sandy soils).

Multivariate analysis on gene expression and environmental factors

700

800

A principal component analysis (PCA) was performed on the gene expression data to investigate factors causing variance in the dataset. The first two principal components (PC) accounted for 68% of the variation in the dataset (44% for the first and 24% for the second PC). The first PC explains the variation between the

two soil types while the second PCA accounts for other variation, including variation induced by differences in land-use. In order to link gene expression differences to the soil environmental factors to which the collembolans were exposed, we performed a CCA. The obvious effect due to soil-type was removed before the analysis. The soil characteristics that showed the highest correlation with gene expression differences were three measures of phosphate: phosphate in pore water (Pw), bio-available phosphate (P-AL) and total phosphate (P-total). Even though these are all phosphate measurements, they showed no correlation between each other, reinforcing that there is a real effect of soil fertility on gene expression regardless of soil-type or land-use. The weight distribution of the genes in the CCA was calculated to test which genes correlated best with increasing phosphate concentration in the soil, either by upregulation or downregulation. Gene correlation weights ranged from )0.053 to 0.050. A cut-off setting of )0.025 and 0.025 was used to exclude genes with a low correlation score, yielding 201 genes that were upregulated with increasing phosphate concentration and 200 genes that were downregulated with increasing phosphate. Both sets of genes were subjected to GO analysis focusing on biological process and molecular function. In the GO terms that were upregulated with increasing phosphate concentrations no clear pattern could be observed. However, in the terms that were negatively correlated with phosphate concentration we found a clear signature of protein synthesis. The terms; ‘translation’ (GO:0006412), ‘transcription from polymerase I, II and III promoters’ (GO:0006360, GO:0006383 and GO:0006367) and ‘RNA elongation from polymerase II promoter’ (GO:0006368) were all significant in biological process.

200

300

Number of juveniles 400 500 600

Discussion

5

10 15 20 Soil arsenic concentration (mg/kg)

25

Fig. 4 The correlation between reproduction (y-axis) and soil arsenic concentration (x-axis). Reproduction in F. candida is negatively affected at higher arsenic concentrations in the soil.  2011 Blackwell Publishing Ltd

In this study we presented a first step towards defining a NOR, Odum et al. 1979) at the level of transcriptional regulation. We investigated the transcriptional plasticity of a soil-dwelling collembolan exposed to different natural soils. The largest source of variation in gene expression patterns was observed between two soil types: sand and clay. It indicates that this collembolan experiences the physicochemical properties of clay and sandy soils as very different from each other, despite the fact that it can survive and reproduce equally well in these soils. F. candida prefers to lay its eggs in the soil rather than on the surface to protect them against predators (Fountain & Hopkin 2005). The texture of the soil is therefore a factor of immediate relevance for egg laying. Additionally, soil texture influences chemical

1152 T . E . D E B O E R E T A L . composition. Clay soils tend to have higher concentrations of metal ions bound to the exchange complex of lutum particles, while sandy soils have lower lutum and organic matter contents with lower concentrations of exchangeable metals. In our soil chemical analysis lutum content was positively correlated with the metals Cr, Zn, Ni and Cu. Despite the fact that it can easily survive in sandy soils and is often found there, F. candida prefers soils that are rich in organic matter (Potapow 2001). In the gene expression GO analysis, we found for sandy soils that GO terms involved in the cell cycle and DNA replication and repair were upregulated while for animals exposed to clay soils GO terms in protein metabolism were upregulated. Collembolans exposed to sandy soils seem to experience more general stress probably due to the lower organic matter content of the soil. However, specific stress response pathways such as those dealing with the metabolism of xenobiotic substances were not found to be regulated in either soil-type. The land-use of the soil had much less impact on gene expression than the soil-type. Still, even though locations were classified to have the same land-use there were differences in, for example, vegetation and fertilization schedules. Of all the land uses, only forest soils induced a specific pattern. Gene transcripts found differentially regulated in forest soils included ABC transporters and glucuronosyl transferases. These genes are involved in the removal of plant secondary metabolites and recalcitrant aromatics which can be harmful to the collembolan. Forest soils are likely to have a higher concentration of humic acids and polyphenols than agricultural soils, due to the higher plant and tree coverage. Decomposition of plant and wood material, such as lignins, by bacteria and fungi can produce residues that are toxic to collembolans when ingested (Berg et al. 2004). The forest soils were the only soils that induced this specific detoxification pathway in F. candida. We also investigated the gene expression data for an interaction effect between soil-type and land-use but no such interaction was found. One of the reasons for this might have been the small number of genes that were differentially expressed between the different land uses. Thus, soil factors associated with land-use do not influence the animal’s physiology to a great extent. Several authors have emphasized the biological relevance of natural variation in gene expression (Whitehead & Crawford 2006). For example, Crawford and Oleksiak (2007) measured expression of metabolic genes in different individuals of an out-bred strain of Fundulus heteroclitus and found that up to 81% of the variation in physiological metabolism in the heart could be explained by gene expression. Heritable variation in gene expression among individuals can also contribute

to genetic differences between populations and evolutionary change (Roelofs et al. 2009). Genetic adaptation might be less dynamic in a parthenogenetic reproducing species such as F. candida. We speculate that plasticity may play a more important role in ecological interactions of this species. However, accumulation of mutations goes faster due to the fact that they are less frequently filtered out by recombination. The 28-days reproduction test showed considerable variation in the number of juveniles among the different soils. F. candida exposed to clay soils showed an average lower (not significant) reproduction but here was no clear pattern of the number of juveniles concerning land-use. However, a significant correlation was established between the concentration of soil arsenic and F. candida fitness, despite the fact that the arsenic concentration in all soils was below the reference value for the Netherlands (29 mg ⁄ kg soil) (VROM 2000). These unexpected sub-lethal effects of arsenic might be due to toxicity of arsenic itself, but it could also be due to an interaction with phosphorus, since arsenic and phosphorus are known to compete during soil sorption and uptake by plants. The high level of variability within the reproduction test may have been caused by the use of unsterilized soils, thus harboUring different microbial communities. We deliberately did not sterilize the soils to keep the soil microbial community intact and to maintain soil ecosystem functions. Kaneda and Kaneko (2002) showed that F. candida growth rate is influenced by soil bacterial activity, which may lead to differential reproductive output after 28 days of exposure. Also small invertebrates such as nematodes, mites and other collembolans were not removed. These biotic factors may have additional effects on the reproduction test over 28 days, but are considered to be less relevant in the 2-day gene expression test. Endogenous F. candida was not observed in the soils. Three natural soils showed less than 50% of F. candida reproduction when compared to reproduction in the standard soil LUFA 2.2 (soils 6, 12 and 14), despite the fact that these soils were all unpolluted. This suggests that soil characteristics alone can confound bioassays and cause effects on F. candida gene expression when the soil is compared to a standard lab soil. Such effects can also occur when testing a suspected polluted soil so that reproduction as an endpoint for soil quality assessment is of limited value and should be supported by additional tests based on gene expression. In a CCA we linked gene expression to the data from the soil chemical analysis. The main soil characteristic correlated with gene expression was soil fertility. Soil fertility is mainly determined by phosphate, which was measured with four different methods (total, soluble,  2011 Blackwell Publishing Ltd

S O I L A R T H R O P O D T R A N S C R I P T I O N A L P L A S T I C I T Y 1153 etc.). Three of the four different methods resulted in strong correlations of phosphate with gene expression. F. candida is often found in agricultural fields and prefers soils with elevated organic matter content (Potapow 2001). The fact that gene expression, even on the short term, correlates with soil fertility indicates that habitat preference of this animal is measurable on the transcriptional level. Soil fertility did not, however, correlate with F. candida reproduction. Higher concentrations of phosphate in the soil were negatively correlated with protein synthesis which indicates that collembolans exposed to lower soil fertility upregulate their protein synthesis. Since the exposure was only 2 days, this means that F. candida, exposed to a (less preferred) soil of low fertility, stages a gene expression cascade to reach homeostasis, which involves increased protein synthesis. In this article we took the first steps to establish a transcriptional NOR for the ecological relevant test animal F. candida by determining transcriptional plasticity under various natural conditions. The soils used in this study have soil types and land uses in which this collembolan is normally found. The F. candida transcriptome seems to be rather plastic, which is indicative of the fact that it can tolerate a wide range of abiotic factors (e.g. soil pH, moisture content, etc.). Plasticity seems to be a common property of soil-living detritivorous and plastic species may be favoured when communities shift under environmental change (Berg & Ellers 2010). We have shown that plasticity comes with significant regulation of gene expression, which may involve 18.5% of the transcriptome. The next step in completing a NOR for this invertebrate would be to vary exposure times and use multiple genotypes from natural populations to investigate response differences over time and in lineages sampled from different soil conditions.

Acknowledgements The authors would like to thank Bart Pieterse and Eiko Kuramae for their help with soil sampling; Francois Rustenburg and Paul Eijk for assisting with the microarray experiments; Ben Nota and Thierry Janssens for help with gene expression analysis and Janine Marie¨n for assistance during RNA isolation. This work was supported by a grant from the Netherlands Genomics Initiative (NGI) to the Ecogenomics consortium ‘Assessing the Living Soil’.

References Alexa A, Rahnenfuhrer J, Lengauer T (2006) Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics, 22, 1600– 1607.

 2011 Blackwell Publishing Ltd

Altman NS, Hua J (2006) Extending the loop design for twochannel microarray experiments. Genetical Research, 88, 153– 163. Bayley M, Holmstrup M (1999) Water vapor absorption in arthropods by accumulation of myoinositol and glucose. Science, 285, 1909–1911. Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate – a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B-Methodological, 57, 289–300. Berg M, Ellers J (2010) Trait plasticity in species interactions: a driving force of community dynamics. Evolutionary Ecology, 24, 617–629. Berg MP, Stoffer M, van den Heuvel HH (2004) Feeding guilds in Collembola based on digestive enzymes. Pedobiologia, 48, 589–601. Cavalieri D, Townsend JP, Hartl DL (2000) Manifold anomalies in gene expression in a vineyard isolate of Saccharomyces cerevisiae revealed by DNA microarray analysis. Proceedings of the National Academy of Sciences of the United States of America, 97, 12369–12374. Crawford DL, Oleksiak MF (2007) The biological importance of measuring individual variation. Journal of Experimental Biology, 210, 1613–1621. Crouau Y, Cazes L (2003) What causes variability in the Folsomia candida reproduction test? Applied Soil Ecology, 22, 175–180. Crouau Y, Moı¨a C (2006) The relative sensitivity of growth and reproduction in the springtail, Folsomia candida, exposed to xenobiotics in the laboratory: an indicator of soil toxicity. Ecotoxicology and Environmental Safety, 64, 115–121. Edgar R, Domrachev M, Lash AE (2002) Gene expression omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Research, 30, 207–210. Edwards D (2003) Non-linear normalization and background correction in one-channel cDNA microarray studies. Bioinformatics, 19, 825–833. Fountain MT, Hopkin SP (2005) Folsomia candida (collembola): A ‘‘Standard’’ soil arthropod. Annual Review of Entomology, 50, 201–222. Gibson G (2002) Microarrays in ecology and evolution: a preview. Molecular Ecology, 11, 17–24. Hattenschwiler S, Tiunov AV, Scheu S (2005) Biodiversity and litter decomposition interrestrial ecosystems. Annual Review of Ecolog,y Evolution and Systematics, 36, 191–218. Huang X, Suyama K, Buchanan J, Zhu AJ, Scott MP (2005) A Drosophila model of the Niemann-Pick type C lysosome storage disease: dnpc1a is required for molting and sterol homeostasis. Development, 132, 5115–5124. ISO (1999) ISO, Soil quality. Inhibition of reproduction of collembola (Folsomia candida). ISO guideline 11267. In: International Standardization Organization, p. 16, Zwitserland. Kammenga JE, Herman MA, Ouborg NJ, Johnson L, Breitling R (2007) Microarray challenges in ecology. Trends in Ecology and Evolution, 22, 273–279. Kaneda S, Kaneko N (2002) Influence of soil quality on the growth of Folsomia candida (Willem) (Collembola). Pedobiologia, 46, 428–439. Kaneda S, Kaneko N (2008) Collembolans feeding on soil affect carbon and nitrogen mineralization by their influence on

1154 T . E . D E B O E R E T A L . microbial and nematode activities. Biology and Fertility of Soils, 44, 435–442. Kerr MK, Martin M, Churchill GA (2000) Analysis of variance for gene expression microarray data. Journal of Computational Biology, 7, 819–837. Noe¨l HL, Hopkin SP, Hutchinson TH, Williams TD, Sibly RM (2006) Population growth rate and carrying capacity for springtails Folsomia Candida exposed to ivermectin. Ecological Applications, 16, 656–665. Nota B, Bosse M, Ylstra B, van Straalen NM, Roelofs D (2009) Transcriptomics reveals extensive inducible biotransformation in the soil-dwelling invertebrate Folsomia candida exposed to phenanthrene. BMC Genomics, 10, 236. Nota B, Timmermans M, Franken O, et al. (2008) Gene expression analysis of collembola in cadmium containing soil. Environmental Science and Technology, 42, 8152–8157. Odum EP, Finn JT, Franz EH (1979) Perturbation theory and the subsidy-stress gradient. Bioscience, 29, 349–352. Oleksiak MF, Churchill GA, Crawford DL (2002) Variation in gene expression within and among natural populations. Nature Genetics, 32, 261–266. Potapow M (2001) Synopses on Palaearctic Collembola, Volume 3. Staatsliches Museum fur Naturkunde Gorlitz, Gorlitz. Pritchard CC, Hsu L, Delrow J, Nelson PS (2001) Project normal: defining normal variance in mouse gene expression. PNAS, 98, 13266–13271. Roelofs D, Aarts M, Schat H, Van Straalen N (2008) Functional ecological genomics to demonstrate general and specific responses to abiotic stress. Functional Ecology, 22, 8–18. Roelofs D, Janssens TKS, Timmermans M, et al. (2009) Adaptive differences in gene expression associated with heavy metal tolerance in the soil arthropod Orchesella cincta. Molecular Ecology, 18, 3227–3239. Saeed AI, Bhagabati NK, Braisted JC, et al. (2006) [9] TM4 microarray software Suite. In: Methods in Enzymology (eds Alan K, Brian O), pp. 134–193. Academic Press, San Diego. Schena M, Shalon D, Davis RW, Brown PO (1995) Quantitative monitoring of gene expression patterns with a complementary DNA microarray. Science, 270, 467–470. Smit CE, Van Gestel CAM (1998) Effects of soil type, prepercolation, and ageing on bioaccumulation and toxicity of zinc for the springtail Folsomia candida. Environmental Toxicology and Chemistry, 17, 1132–1141. Smyth GK (2004) Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology, 3, Issue 1, Article 3. Timmermans M, de Boer M, Nota B, et al. (2007) Collembase: a repository for springtail genomics and soil quality assessment. BMC Genomics, 8, 341. Timmermans MJTN, Roelofs D, Nota B, Ylstra B, Holmstrup M (2009) Sugar sweet springtails: on the transcriptional

response of Folsomia candida (Collembola) to desiccation stress. Insect Molecular Biology, 18, 737–746. Van Straalen NM, Roelofs D (2008) Genomics technology for assessing soil pollution. Journal of Biology, 7, 19. VROM (2000) Streefwaarden en interventiewaarden bodemsanering. Staatscourant, 39, 8–9. Whitehead A, Crawford DL (2006) Variation within and among species in gene expression: raw material for evolution. Molecular Ecology, 15, 1197–1211. Witten DM, Tibshirani R, Hastie T (2009) A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis. Biostatistics, 10, 515–534.

We are interested in the molecular mechanisms underlying stress responses in soil invertebrates. These molecular responses to stress may be used as biomarkers to measure the impact of chemicals and pollutants. We mainly focus on gene expression analysis of collembolans that are exposed to chemically spiked soils and polluted field soils. To be able to distinguish between genomic responses induced by pollution and natural gene expression responses, we are developing a reference database that consists of gene expression profiles of collembolans exposed to unpolluted field soils.

Supporting information Additional supporting information may be found in the online version of this article. Table S1 List of significant genes, description, fold change and adjusted P-value encountered in the MAANOVA analysis between the different soil types, clay and sand. Descriptions were derived from BLAST analysis with a minimal e-value of 10)10. Positive fold change means upregulated in clay soils. Fig. S1 The microarray design that was used in this study. The design consisted of two interwoven loops. One loop was used for the clay soil samples and one loop for the sandy soil samples. The circles represent the soils while the arrows represent the two colour microarrays (the point of the arrow represents the Cy5 sample while the tail represents the Cy3 sample). Four replicates were used per soil sample; two labelled with Cy3 and two with Cy5. No replicates of the same soils were every hybridized against the same soil more than once. LCz and LCc stand for the LUFA soil controls incorporated in the design. Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

 2011 Blackwell Publishing Ltd