Seminal fluid proteins (Sfps) alter female behaviour and physiology and can mediate sexual conflict. In Drosophila melanogaster, a single Sfp, the sex peptide (SP), triggers remarkable post-mating responses in females, including altered fecundity, feeding, immunity and sexual receptivity. These effects can favour the evolutionary interests of males while generating costs in females. We tested the hypothesis that SP is an upstream master-regulator able to induce diverse phenotypes through efficient induction of widespread transcriptional changes in females. We profiled mRNA responses to SP in adult female abdomen (Abd) and head+thorax (HT) tissues using microarrays at 3 and 6 h following mating. SP elicited a rich, subtle signature of temporally and spatially controlled mRNAs. There were significant alterations to genes linked to egg development, early embryogenesis, immunity, nutrient sensing, behaviour and, unexpectedly, phototransduction. There was substantially more variation in the direction of differential expression across time points in the HT versus Abd. The results support the idea that SP is an important regulator of gene expression in females. The expression of many genes in one sex can therefore be under the influence of a regulator expressed in the other. This could influence the extent of sexual conflict both within and between loci.
The interplay between the sexes over reproduction generates substantial phenotypic variability upon which selection can act [1–3]. To understand the pace, dynamics and trajectory of coevolution arising from this interplay, it is necessary to understand the molecular interactions between males and females . Costs and benefits of reproductive interactions will determine whether coevolution is antagonistic [5–7]. A mechanistic understanding of costs should prove particularly informative in revealing the nature and extent of the potential manipulation of gene expression in one sex by the other. This would allow an assessment of the potential ease with which females could evolve resistance to manipulative adaptations  and reveal any constraints arising from the genetic architecture of the pathways involved.
Insights into the molecular interplay between males and females have been gained through economic studies of single genes [2,8,9], studies of fitness effects of interactions between different male and female genotypes [1,10], the characterization of male and female reproductive proteins [3,11] and the measurement of genome-wide responses to mating [12–17]. Of these methods, only genome-wide studies can identify the overall extent to which male and female transcriptomes are altered when the sexes interact, and they offer a valuable opportunity to probe the underlying basis of the interacting pathways involved. Much of the mechanistic work has used the powerful genetic resources of the fruitfly Drosophila melanogaster. Particularly well studied are the functions and fitness effects of seminal fluid proteins (Sfps) in males [2,18], which cause a substantial reprogramming of female behaviour and physiology . The genome-wide profiles of gene expression changes in response to courtship and mating [12,13,16,17] and to the receipt of Sfps [12,15,19] have been investigated. Here, we aimed to capture an in-depth genomic profile of gene expression changes resulting from the transfer of a single Sfp with striking fitness effects [20–26], the sex peptide (SP).
Sex peptide can produce a startling array of female responses, e.g. reduced female receptivity, increased egg laying, feeding and antimicrobial peptide production, and altered sleep patterns and water balance [20–26]. It also adheres to sperm and influences the rate of sperm release from storage [27,28]. The receptor for SP  is expressed in the female brain, nervous system and reproductive system, and responses to SP depend upon a circuit of neurons that express the sex determination genes doublesex  and fruitless . Insight into the extent of the influence of SP in females is of interest because of the finding that SP can mediate sexual conflict. Females exposed to high levels of SP have reduced fitness , the magnitude of which is dependent upon female nutritional status . Males, on the other hand, benefit from SP transfer , an effect that is also nutrition dependent . These findings show that costs arising from sexual conflict are context dependent, and highlight the role of nutritional status in determining their magnitude. They are also consistent with the hypothesis that males can manipulate female reproductive processes by altering the sensitivity of nutrient-sensing pathways. Our aim in this study was to reveal the full extent of genomic responses of female D. melanogaster to SP, to measure at the genomic level the extent of male manipulation and to predict potential evolutionary genomic responses in females.
To examine differences in mRNA expression in females following SP receipt, we mated wild-type females to males lacking SP (null, SP0) or to control, SP-producing (SP+), males . We assayed mRNAs regulated by SP at 3 and 6 h after mating, to coincide with phenotypic responses observed at 3 h [20,21,23,33–35] and through 6 h and beyond [14,20]. We conducted fourfold independent biological replication. We profiled abdomen (Abd) and head and thorax (HT) tissues separately. This was to minimize the loss of mRNA signal arising from tissue-specific expression  and to capture distinct signals arising from the reproductive tissue (Abd) versus central nervous system (HT). Our goal was to identify the full range of biological processes affected by SP receipt. Hence, we opted for a gene category enrichment analysis (Catmap ), rather than defining lists of differentially expressed genes with relatively low fold changes based on a cutoff. All expressed genes were ranked based on significance (Bay-t statistic, see below), and this list was then analysed for over-represented functional annotation terms. This Catmap approach integrated information from all expressed genes, giving sensitive detection of, and well-supported biological insight into, process-level changes induced by receipt of SP.
(b) Fly stocks and husbandry
Females were from an outbred, laboratory-adapted Dahomey stock maintained in population cages with overlapping generations. Control and SP null males were derived as in . SP0 bears a non-functional SP gene opposite Δ130, a deletion that covers SP . The SP+ genetically matched control gives wild-type SP expression , but contains both the mutant and wild-type SP genes. We backcrossed the deletion and mutant alleles into the Dahomey genetic background (as in ). Flies were kept at 25°C in non-humidified rooms on a 12 L : 12 D cycle and supplied throughout with sugar yeast (SY) food (as in ) supplemented with excess live yeast granules.
(c) Matings with SP0 and SP+ males
Parents of experimental flies were 4 days old when their eggs were sampled. All flies were grown at a low standard density (18 µl eggs/70 ml food ). Virgin Dahomey females were collected on ice within 6 h of eclosion and stored in single sex groups of 10. SP0 and SP+ males were collected within 24 h of eclosion using ice or CO2 anaesthesia, and stored in single sex groups of 10. Mating assays were performed when females were 6 and males were 5–6 days old. The day before matings, two males were aspirated into each mating vial (Nunc cryotube, 4.5 ml with starfoot, containing SY food and a standard-size, live yeast grain). On the day of the assay, a single female was aspirated into each vial (the first female into an SP0 vial, the second into a SP+, the third into an SP0 vial and so on). Females that mated for less than 10 min were discarded, as complete ejaculate transfer was unlikely . Vials containing once-mated females were formed into matched pairs, in which SP0 and SP+ matings began within 5 min of each other. Following matings, males were removed. Profiling time points were 3 h and 6 h after the start of mating (ASM). The first pair of females was allocated to the 3 h profile, the second to the 6 h, the third to the 3 h and so on. At 3 and 6 h ASM, the females were frozen in their cryotubes in liquid N2. After freezing but before RNA extraction, Abd tissues were separated from HT without defrosting. Body parts (HT or Abd) of females within each SP treatment (SP0 or SP+) within each time point were combined into a single microcentrifuge tube (eight tubes in total). The experiment was replicated, using exactly the same procedure, four times. Thus, a total of 32 tubes were produced (i.e. two tissue types (HT, Abd) × 2 treatments (SP+, SP0) × 2 time points (3 h, 6 h) × 4 replicates = 32). The final sample size of body parts within each tube was N = 61–74. To verify the phenotype of backcrossed SP0 and SP+ lines, we performed a bioassay. During replicate 2, we performed 20 additional matings for each treatment. Females were transferred individually to fresh vials after mating and allowed to oviposit for 24 h. They were then introduced to vials containing two wild-type males, and the number of females that remated within 30 min was recorded. As expected, females mated to SP0 males laid significantly fewer eggs than controls (F1,38 = 14.28, p = 0.0006) and were significantly more receptive (, p < 0.0001, data not shown).
(d) RNA extractions, labelling, hybridizations
Total RNA (TRNA) was isolated using the mirVana kit (Ambion) following the manufacturer's protocol. Quality and quantity of the isolated RNA were verified using spectrophotometry and gel electrophoresis. The TRNA sample was used for mRNA expression analysis and the remainder stored at −80°C for qRT-PCR validation. mRNA array assays were conducted at the Glasgow Polyomics Facility (University of Glasgow). Biotinylated cRNA probes were generated with the one-cycle Genechip target labelling protocol and hybridized to high-density oligonucleotide Affymetrix microarrays (Drosophila Genome 2.0 Array), as described in the Affymetrix GeneChip Expression Analysis Technical Manual (2000).
(e) Gene expression (mRNA) data analyses
Microarray data are lodged in ArrayExpress (www.ebi.ac.uk/arrayexpress; accession nos. E-MEXP-2070 and E-MEXP-2074 for the HT and Abd datasets, respectively). With the exception of Catmap perl program for gene category analysis, all analyses were done using R v. 2.6.2  and BioConductor . Preliminary quality control identified one HT outlier microarray, which was excluded from further analyses. Hence, a total of 16 arrays were analysed for the Abd and 15 for the HT. Hierarchical clustering among samples and ANOVA analysis indicated a batching effect correlated with the ‘biological replicate’ random variable. ANOVA revealed that variance among time-points was heterogeneous, suggesting that time-points should be treated separately.
The numbers of present/expressed genes in the HT and Abd datasets were determined using MAS5 present/absent calls, based on GC-RMA-transformed PM and MM values, as previously implemented in R . A transcript was considered present if the p-value was <0.111 at either time-point/genotype. Normalizations were performed using GC-RMA on all probes , followed by Loess  on present probe sets. We included the biological replication batching structure in global statistical analysis by conducting paired comparisons. The Cyber-T  Bayesian statistic was used to perform paired t-tests between SP+ and SP0 strains for each time-point. False discovery rates were estimated using the Storey method [46–48], as implemented in R (Qvalue package).
To identify significant differential expression of functionally related categories of genes, we used Catmap analysis . This program was populated with functional annotations for the genome of D. melanogaster, using universal vocabularies from GeneOntology (GO) 1.107 , Interpro  and KEGG v. 44 [51–53], as well as links to PubMed identifiers for published genes retrieved from Flybase . An additional functional category comprised groups of genes with known tissue-specific expression . This allowed us to test for overlap of particular categories of genes with clusters of genes expressed in specific tissues—hence to discover with greater resolution in which tissues the Abd and HT genes were differentially expressed. Raw data (.CEL) files for 15 fly tissues from FlyAtlas  (GEO accession no. GSE7763) were used together with an additional dataset (N. Alic, M. E. Giannakou & L. Partridge 2009, unpublished data) corresponding to adult cuticle and associated fat body in white Dahomey S1106/+ flies . The number of present genes in these datasets was determined as described above , and normalizations performed on present genes using both GC-RMA  with a Loess  option, and MAS v. 5.0  with a quantiles  option. Soft clustering of these FlyAtlas-expressed genes was achieved using the mfuzz program [59,60] to incorporate the resulting clusters.
For each time-point and body part, Catmap analysis was conducted on two lists of genes ranked based on the Bayes t statistic from Cyber-T. One was based on ranking all expressed genes from the most significantly up- to downregulated, and the second on the reverse. As a category input file, the same list of present genes mapped to functional categories was supplied. As Catmap p-values are positively biased towards categories with many genes under the gene permutation null hypothesis , categories comprising less than five genes were discarded. The one-sided Wilcoxon rank sum was used to generate a score based on the sum of the rankings of present genes with a particular functional annotation. The significance of that score (the p-value) was calculated analytically based on a random gene rank distribution. Adjusted p-values (correction for multiple testing) were calculated based on 1000 random permutations of gene rankings. Following identification of enriched functional categories by Catmap, we tested the significance of functional overlaps between tissue-specific FlyAtlas clusters and other functional categories (such as GO, KEGG, etc.) that showed similar regulation profiles, using Fisher exact tests.
(g) Quantitative RT-PCR
To confirm the up- and downregulation of gene categories, individual genes were selected for qRT-PCR. We chose 10 genes from 3 h in the HT, and 10 from 6 h in the Abd. We used the same samples for qRT-PCR validation as for the microarrays. The choice was based on the degree of up- or downregulation (expression fold change) and statistical confidence (ranked p-values from the Bayes t statistic  and the q-value after Storey multiple testing correction [47,48]). RNA was quantified using a NanoDrop ND1000 spectrophotometer (LabTech, Ringmer UK). One microgram of TRNA was reverse-transcribed using 200 ng Random Primers and Supercript ll reverse transcriptase (Invitrogen, Paisley, UK), according to the manufacturer's instructions. qRT-PCR reactions were performed using the Taqman ABI 7900, with 12.5 µl of SYBR Green PCR Master Mix (Applied Biosystems, Warrington), 200 nM of each primer and 5 ng of cDNA in 25 µl total volume. Tests with control genes gave extremely low variability, and technical replicates were not used. Dissociation curves were run following each reaction to ensure single product amplification. No contamination in the melting curves was found for any primers used (for sequences, see the electronic supplementary material, table S1). To determine relative RNA levels, standard curve analysis was used. qRT-PCR reactions were prepared by using the cDNA from one sample and making 2-fold serial dilutions covering the range equivalent to 100 ng to 3.125 ng of RNA. αtubulin84β (FBgn0003884) was used as an endogeneous control, as this was more stable than 18S rRNA.
3. Results and discussion
Receipt of SP led to widespread but subtle changes in mRNAs (see figure 1 and the electronic supplementary material, figures S1–S9 and table S2). Throughout, upregulated refers to functional categories (groups) of genes with higher expression following receipt of SP, and downregulated to lowered expression induced by SP, in comparison with SP-lacking controls.
(a) mRNA responses: global patterns
Changes in gene expression were globally small in expression ratios (see the electronic supplementary material, figures S8, S9 and table S3), consistent with previous studies [12,13]. We assessed both quantitative (correlation) and qualitative (direction) validation of the microarray data by qRT-PCR. This was consistent overall with the differential expression detected by microarrays (see the electronic supplementary material, figures S8, S9 and table S3). Quantitative validation was provided by the significant positive correlation (r = 0.668, p = 0.002, n = 19) between fold changes from Affymetrix arrays versus qRT-PCR (see the electronic supplementary material, figure S9). Of the 20 genes tested by qRT-PCR, 18 showed directional agreement with the array data (see the electronic supplementary material, table S3 and figure S8; binomial probability p = 0.0004). However, the small fold changes observed render statistical (e.g. t-test) analysis of the qRT-PCR data insensitive (particularly for fold changes less than 1.4 ). Consistent with this, few genes were confirmed by t-tests as significantly regulated (see the electronic supplementary material, table S3 and figure S8). Exact regulation patterns for genes subjected to qRT-PCR should therefore be treated with caution. Indeed, this was a motivation for exploring expression patterns through the process-level Catmap analysis, instead of focusing on individual genes drawn from gene lists.
A summary of Catmap analysis of all expressed genes in the HT (8274 genes) and Abd (7626 genes) (figure 1; full results in the electronic supplementary material, table S2) illustrates enriched categories from gene functional annotations . In the Abd, we observed subtle but widespread downregulation of mRNAs at 3 h, a pattern that persisted but was less marked at 6 h. Functional enrichment analysis illustrated that 452 and 199 gene categories were downregulated at 3 and 6 h, respectively (see figure 1 and electronic supplementary material, figures S1–S3). Only a few categories showed significant upregulation, i.e. 68 and 26 categories at 3 and 6 h. In contrast, the HT showed a much more dynamic pattern. At 3 h, 97 categories were up- and 134 downregulated; at 6 h, 73 categories were up- and only 19 downregulated. It was common in the HT profile to observe significant enrichment within a category at only one time point or even opposite profiles between 3 and 6 h (see figure 1 and the electronic supplementary material, figures S4–S6).
mRNA responses to SP in the HT encompassed diverse biological processes, and the direction of differential expression varied markedly within and across timepoints (figure 1a). In contrast, Abd profiles exhibited a more consistent pattern of downregulation, with changes in genes related to a narrower set of potential functions. There were several specific examples of tight regulation across time. For example, differential expression of genes regulated by TOR (target of rapamycin), a nutrient-sensitive regulator of cell growth, proliferation, dietary choices and lifespan , was significantly downregulated at 3 h in the HT, but showed an opposite expression profile at 6 h (see the electronic supplementary material, figure S5). This finding suggests precise transcriptional control over short time-scales. There were also examples of apparent systemic coordination of classes of genes across different body parts. For example, Jonah genes and peptidases showed time- and tissue-specific activation and were significantly upregulated at 3 h in the Abd and downregulated in the HT at 6 h (see figure 1 and the electronic supplementary material, figures S2, S5). We hypothesize that these proteins could have a role in digestion, or play a role in defense in the midgut.
(b) mRNA responses: potential functions of differentially expressed mRNAs
In the HT, there was significant differential expression in gene categories related to signalling pathways, behaviour, cell communication and development (figure 1). Novel and unanticipated signatures included the upregulation of genes involved in phototransduction , highly time-specific regulation of genes regulated by TOR and downregulation of chorion genes (see the electronic supplementary material, figure S5) . In marked contrast to the HT, Abd profiles showed changes in genes related to fewer potential functions. We describe major mRNA signatures below.
In the HT profile at 3 h post-mating, there was significant upregulation of processes related to phototransduction (figure 1a). Related components, such as the G-signalling pathway, transmembrane receptor and molecular transducer activities were all differentially expressed, as well as genes involved in cell communication, brain and behaviour, and synaptic transmissions (see figure 1 and the electronic supplementary material, figures S4–S6). Lysosomal tetraspanin, a key protein related to phototransduction, was also significantly upregulated at 3 h (see the electronic supplementary material, figure S5). To determine the specific tissues in which these sets of genes were expressed, we cross-referenced to gene clusters built from the FlyAtlas data . Two clusters upregulated at 3 h (11 and 19, comprising genes expressed in the brain and in the thoracic ganglion, respectively), contained genes with significant functional overlap with genes upregulated in response to the SP at 3 h (see figure 1a and the electronic supplementary material, figure S7a; Fisher test p-value < 10–7). The overlapping functional enrichment was for phototransduction, G-signalling, neuropeptide function and synapse structure (cluster 19), as well as for neurological processes, cell surface and molecular transducers (cluster 11). This result suggests that visual perception is altered following the receipt of SP, perhaps to facilitate the search for suitable oviposition sites.
(ii) TOR pathway
The amino acid sensing TOR pathway was downregulated by SP in the HT at 3 h post mating, along with related categories in energy metabolism and energy biogenesis (see the electronic supplementary material, figure S5). This suggests a link between the receipt of SP and the regulation of nutrient sensing [66,67], which could mediate the known role of SP in feeding and nutrient balancing [24,68]. This is consistent with the finding that the effects of Sfps, and the costs of SP receipt, are significantly altered by dietary manipulations [6,32]. Notably, the TOR pathway showed the opposite profile (upregulation) at 6 h, indicating that SP induces rapid changes in nutrient sensing. Signalling related to increased feeding to support egg production is potentially induced relatively late (from 6 h).
(iii) Ovary and development genes
The majority of categories differentially expressed in the Abd were involved in diverse aspects of egg and early embryo development (see figure 1 and the electronic supplementary material, figure S2). Examples include cell motility and migration, cytoskeleton genes, signal transduction (e.g. JAK/STAT signalling), regulation of transcription and RNA localization, regulation of divisions/cell cycle and cell size and gamete formation (see figure 1 and the electronic supplementary material, figures S1–S3). These signals reflect the role of SP in activating egg development (e.g. cytoskeleton genes) and also the increased rate of production of eggs loaded with translationally inactive maternal mRNAs that direct early embryonic development [69,70]. Early development genes could therefore be differentially expressed in females because of transcription of maternally loaded mRNAs, and/or from early zygotic transcription in retained fertilized eggs within the female reproductive tract. However, the majority of functional categories were significantly down- rather than upregulated in the Abd, in contrast to the expectation given that receipt of SP significantly elevates egg production [20,21]. A plausible explanation is the loss of maternal mRNAs [69,70] in the significantly increased numbers of eggs laid. An additional and unexpected signature was downregulation in the HT of ‘ovary-related’ (e.g. chorion) genes (see figure 1 and the electronic supplementary material, figures S4–S6). These did not show significant overlap with the Abd profiles noted earlier. Significant downregulation was observed for the maternal gene oskar, various cyclins, and genes involved in mRNA translation, splicing and nucleic acid binding via Myc-Max-Mad/Mnt transcription factors (see figure 1 and the electronic supplementary material, figure S5). Such changes related to the control of growth and proliferation are predicted for egg development and maturation , but are expected for Abd and not for HT profiles. By cross-referencing to clusters of tissue-specific expression, we found a significant overlap with a set of genes expressed in the cuticle-fat body as well as in the ovaries (cluster 15, figure 1a; electronic supplementary material, figure S7a). The enrichment of these genes in the fat body potentially resolves the paradox of ‘ovary genes’ in the HT (also observed previously by others ). Consistent with this, reduced expression of CG1478 (Chorion protein 36 or Cp36) in the HT was supported by the qRT-PCR (see table 1 and the electronic supplementary material, figure S8a).
(iv) Immune genes
As expected [12,13,19], immune genes formed part of SP responses. They were one of the few categories to be significantly upregulated in the Abd. At 3 h, there was evidence for differential expression in peptidases, proteases/endopeptidases, hydrolases, in Toll-lmd-mediated immune responses and defence-related Jonah genes (e.g. electronic supplementary material, figures S2, S5). There was significant overlap with genes specifically expressed in the midgut (cluster 7, electronic supplementary material, figure S7b) and serine-type proteases, immune response and Jonah genes (Fisher test p ≤ 5×10–10). Part of the immune response was also mediated via a cluster of genes upregulated in the fat body (cluster 15, figure 1b; electronic supplementary material, figure S7b), significantly correlated with Toll-lmd pathway, oxidoreductases and serine proteases (p ≤ 7×10–4). A strong association with Toll-lmd pathway, immune response (p ≤ 3×10–21) and endopeptidase genes (p = 8.8×10–10) was also found in differentially expressed spermathecae genes (cluster 10, figure 1b; electronic supplementary material, figure S7b). This cluster displayed a continuous pattern of upregulation at both 3 and 6 h, suggesting a sustained activation of immune responses in mated females. Novel aspects of the immune response therefore became apparent by the cross-referencing to Flyatlas . First, SP induced a defence response in genes expressed in the gut (e.g. Jonah genes). Second, SP triggered a systemic immune response involving Toll-Imd, as noted previously . Third, evidence that immune genes were also expressed in tissues such as the spermathecae suggests that there was also a local immune response . These data therefore reveal a multi-faceted immune response to SP.
(c) Signatures of sexual conflict
The extent of female transcriptional responses suggests that SP is an upstream master regulator of many post-mating processes. This reveals multiple routes by which females might alter their responses to SP. It also suggests that global, evolutionary responses by females to SP are likely to be complex. Such responses could be constrained through pleiotropic effects on essential processes ; for example, SP signalling depends upon neural circuits that express sex-determination genes . Hence, SP is plugged into fundamental biological processes of the female, and dampening down global SP responses might not be feasible. Opportunities for antagonistic selection depend upon the magnitude of the difference in the optimum value of reproductive traits across the sexes, and mechanisms available in each sex to enforce their interests . It has been suggested that the evolution of sex-limited expression can release genes from intra-locus sexual conflict by removing the tug of war over optimum expression levels . However, such loci can then quickly/simultaneously become influenced by other loci (e.g. such as SP) and become subject instead to inter-locus conflict. Hence, the evolution of sex limitation per se does not represent evidence of conflict resolution in general (whether within or between loci). SP can also influence the expression of loci in females that are both sex-limited and not sex-limited and could, in principle, simultaneously influence the extent of conflict between loci subject to both intra- and inter-locus sexual conflict. Table 1 summarizes expression signatures from this study consistent with traits affected by SP, and predicts opportunities for antagonistic selection. It is possible that traits over which there is conflict could be typified by contrasting expression profiles between virgins and mated females that do and do not receive SP. The transient expression profiles we observed for TOR signalling along the time course of mating (see the electronic supplementary material, figure S5) could represent such an example, but expression profiles of virgins are needed to confirm this. Many changes in response to SP were in genes expressed in both sexes that have contrasting (often essential) roles at different times (table 1). To further explore how the processes in table 1 might reflect sexual conflict at the transcriptional level, one could compare and manipulate the expression of such genes in both sexes at different life-history stages.
Receipt of SP caused significant differential expression in mRNAs. The rich mRNA signatures reflected known and new responses in previously unknown pathways (e.g. phototransduction). Responses in the HT were varied across different gene categories and time. Complexity was evident in several examples of tight temporal and spatial control of mRNA expression. The results show the widespread and complex nature of SP-induced alteration of the female transcriptome. Hence, females that are under selection to alter their responses to the amount and/or potency of SP may have a formidable set of ‘obstacles’ in their way, as would be predicted, given the sexually antagonistic effects of SP [6,18].
We thank Andy Barnes, James Boone, Su Davies, Yasmine Driege, Sebastian Grönke and Giovanna Vinti for assisting with the experiments. We thank Daniela Wieser for help with data analysis. Funding was provided by the NERC, BBSRC, The Royal Society, The Lloyd's Tercentenary Foundation and the Wellcome Trust.
- Received July 18, 2012.
- Accepted August 20, 2012.
- This journal is © 2012 The Royal Society