Hsp90-downregulation influences the heat-shock response, innate immune response and onset of oocyte development in nematodes

Hsp90 is a molecular chaperone involved in the regulation and maturation of kinases and transcription factors. In Caenorhabditis elegans, it contributes to the development of fertility, maintenance of muscle structure, the regulation of heat-shock response and dauer state. To understand the consequences of Hsp90-depletion, we studied Hsp90 RNAi-treated nematodes by DNA microarrays and mass spectrometry. We find that upon development of phenotypes the levels of chaperones and Hsp90 cofactors are increased, while specific proteins related to the innate immune response are depleted. In microarrays, we further find many differentially expressed genes related to gonad and larval development. These genes form an expression cluster that is regulated independently from the immune response implying separate pathways of Hsp90-involvement. Using fluorescent reporter strains for the differentially expressed immune response genes skr-5, dod-24 and clec-60 we observe that their activity in intestinal tissues is influenced by Hsp90-depletion. Instead, effects on the development are evident in both gonad arms. After Hsp90-depletion, changes can be observed in early embryos and adults containing fluorescence-tagged versions of SEPA-1, CAV-1 or PUD-1, all of which are downregulated after Hsp90-depletion. Our observations identify molecular events for Hsp90-RNAi induced phenotypes during development and immune responses, which may help to separately investigate independent Hsp90-influenced processes that are relevant during the nematode’s life and development.


Introduction
Hsp90 is an ATP-dependent molecular chaperone conserved from bacteria to mammals. The cytosolic version of Hsp90 is essential in all eukaryotes. In the nematode C. elegans the Hsp90-homologue DAF-21 is required for development [1][2][3]. This is evident from the complex phenotype that results from RNA interference (RNAi) against Hsp90, which combines

Mass spectrometric analysis
To perform a mass spectrometric analysis the RNAi treatment was performed on plates made from bacterial M9-medium supplemented with ampicillin, tetracycline and 1 mM IPTG. This bacterial growth medium contains 2 g/L glucose and 1 g/L NH 4 SO 4 as the carbon and nitrogen sources [24,25]. To label the nematodes with 15 N, 1 g/L 15 NH 4 SO 4 (euriso-top GmbH, Saarbrücken, Germany) was used instead of the normal 14 NH 4 SO 4 . Nematode growth was not affected by the 15 N-containing nitrogen source. Nematodes were harvested after 2.5 days by washing them off the plates with M9-buffer. They were separated from residual bacteria by gravity sedimentation in a 15 mL plastic tube. The supernatant was removed and after three such washing steps the nematodes were lysed in a Retsch Mixer Mill MM400 in 1 mL M9buffer. The lysate was cleared by centrifugation and the protein concentration was determined by Bradford reagent. The isotope-tagged lysates were then mixed at the calculated ratio to ensure that equal amounts of each sample were combined. This sample was frozen and stored at -80˚C until mass spectrometric analysis was performed. After a tryptic digest of the proteins, 15 μL of the sample peptides were analyzed by nanoLC-MS/MS on an UltiMate1 3000 UPLC system (Thermo Fisher Scientific) online coupled to a Q Exactive Hybrid-Quadropole-Orbitrap mass spectrometer (Thermo Fisher Scientific) with electrospray ion source (ESI). Initially, the samples were desalted and concentrated on a 75 μm x 20 mm Acclaim1 PepMap100 C18 column with 5 μm particle size (Thermo Fisher Scientific) with a flow rate of 30 μL/min using 95% solvent A (0.1% TFA in H 2 O dest ) and 5% solvent B (0.1% TFA, 50% ACN in H 2 O dest ) for 7 minutes. Peptides were subsequently loaded onto a 75 μm x 50 cm Acclaim1 PepMap100 RSLC C18 column with 2 μm particle size (Thermo Fisher Scientific) using a mixture of 95% solvent C (0.1% formic acid (FA) in H 2 O dest ) and 5% solvent D (0.1% FA, 84% ACN in H 2 O dest ). The elution of peptides was performed at a constant flow rate of 400 nL/min in a linear gradient of 5% to 40% solvent D over 120 minutes. Mass spectrometry was performed in the data-dependent acquisition (DDA) mode with Full scan Fourier transform mass spectrometry (FTMS) acquired in an m/z range of 350 to 1400 with a resolution of 70,000. The ten most intense peptides (charge range ! +2) in the FTMS-scan were selected for higher-energy collisional dissociation (HCD) with the collision energy (NCE) set to 27. The dynamic exclusion time was set on 30 seconds with an MS/ MS resolution of 35.000 within an m/z range of 300 to 2000. Data analysis was performed with Mascot Distiller (Matrix Science Inc., USA) as described [25] using a Uniprot database for C. elegans. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD006185 and 10.6019/PXD006185 [PubMed ID: 26527722].

Microarray sample preparation and handling
After RNAi treatment for roughly 2.5 days nematodes were harvested for DNA microarray experiments. Hsp90-depleted and control vector treated nematodes were washed off the plates and collected in 15 mL plastic tubes. Nematodes were separated from residual bacteria by gravity sedimentation. The worms were frozen at -80˚C. RNA extraction and analysis of the mRNA levels was performed at the Kompetenzzentrum für Fluoreszente Bioanalytik (KFB) as a for-fee service. For this analysis Affymetrix GeneChip C. elegans Genome Arrays were used and the data for all genes were RMA-normalized. All microarray datasets will be uploaded to the GEO repository.

Network calculation and evaluation
To compare expression values under various conditions and to identify the most prominent hits we used Microsoft Excel and the clustering program ClusterEx, which was implemented by us to perform semi-quantitative gene clustering based on correlations found in preformed coexpression databases [26,27]. To obtain the initial hit list the differential expression values from the three biological replicates were averaged to yield the mean expression difference and the standard deviations for each gene. The 250 strongest affected genes for each direction were retained and non-protein coding genes were excluded. To construct networks we obtained the 40 highest-ranking coregulated partner genes for each hit from a C. elegans coexpression database obtained from http://seek.princeton.edu/modSeek/worm [28]. Coexpression counts for the hit genes were collected by ClusterEx from this database [26], adding 41x41 connections per included hit gene, leading to network matrices composed of 420250 connections for 250 genes. In network calculations based on real expression data most hits were connected to many other hit genes. This was not the case for random control groups of genes, where few hit-to-hit pairs were generated, but for most genes no coexpressed hit was included. To test the quality of the clusters in the networks, further genes, which had most connections to hit genes, were determined, integrated into the networks (marked with a pink border) and evaluated to confirm the predictive accuracy (CoRegScore) of the network/database setup. This procedure was automated in ClusterEx as described [26]: In short, all genes were ranked according to their expression changes and the positions of the predicted genes were determined and summed up. This value then was used to calculate the CoRegScore between 100 and -100 as described [26]: In the following equation the sum of the predicted genes (AddedRealPositions) is set in relation to the best possible prediction (AddedIdealPositions) and the worst possible prediction (AddedWorstPositions).
For statistical analysis random controls were determined. To this end 100 similar-sized random lists of genes were generated and the analysis was performed in an identical manner. Based on the resulting mean and standard deviation, the significance/p-value for each parameter was determined by calculating the Z-score (http://www.socscistatistics.com/tests/ztest/ zscorecalculator.aspx) and using a one-tailed hypothesis at http://www.socscistatistics.com/ pvalues/normaldistribution.aspx. For all networks in this study the connection parameters show a high level of significance compared to random gene lists (p<0.001).
To test the influence of the specific database, we performed identical calculations with coexpression relationships downloaded from SPELL (http://spell.caltech.edu:3000/) [29] or from COXPRESdb [30] and qualitatively similar networks were obtained. To allow other users to use the described clustering approach for their own C. elegans expression data, we will include it into our webserver at www.clusterex.de [26].

Network visualization and layouting
Networks were finally visualized in CytoScape based on the ClusterEx output files containing all the hit-to-hit connections and their occurrence counts. The final layout was based on the "Edge-weighted Spring Embedded Layout", which clusters the genes according to their connections and brings highly connected genes close to each other [31,32]. Individual nodes were moved only to prevent overlapping nodes in plot areas, where clustering was very high, but care was taken not to alter the location relative to neighboring nodes.
Network construction for the mass spectrometry samples was also performed with Clus-terEx. As no database reporting on coexpression at the protein level was available to us, we also used the modSEEK transcriptional coexpression database here. To compensate for the very small number of hits and the many proteins, for which we did not obtain quantitative expression values in the mass spectrometry experiment, we used the top 100 coregulated genes to construct the networks, instead of 40 as for the microarray data.

Comparison of different microarray experiments regarding similarities
Data from our arrays were compared to publicly available experiments investigating RNAi against sbp-1 and sams-1 (GSM1816551-GSM1816559) [33]. The corresponding CEL-files were obtained from the GEO repository [34]. Expression values were extracted and normalized with RMADataConv and RMAExpress [35]. Similarly we obtained the CEL-files and result files for the response of N2 nematodes to pathogenic bacteria, in particular to the Vibrio cholerae strains VC109 and VC110 (GSE34026) [36] and to Pseudomonas aeruginosa (GSE5 793) [37,38]. To visualize the similarity to the Hsp90 RNAi response, the expression values for the Hsp90-RNAi responsive genes were retrieved from the experiments and then used to color the Hsp90-depletion response according to the respective data. The same color code was used in the supplemental figures to allow a comparison of the same genes within the visualized networks.
To quantify the similarity with our responses we tested, to what extent the hits from Hsp90-responsive clusters are upregulated in these publicly available experiments. To this end we implemented a routine, which calculates a score (UpRegScore) on the scale from 100 (perfect up-regulation) to -100 (perfect down-regulation). It takes as input any list of genes, e.g. from a cluster to be tested, and determines to what extent each gene of this set is upregulated in a chosen microarray experiment. The UpRegScore is calculated based on the positions of these genes in a list of all genes ranked according to their differential expression with the most upregulated on top. It then compares the summed up positions for the investigated genes (AddedRealPositions) to the best possible (AddedIdealUpPositions) and worst possible scenarios (AddedIdealDownPositions), resulting in positive values up to +100 if the group is highly upregulated or negative scores if it is downregulated: AddedRealPositions À AddedIdealUpPositions AddedIdealDownPositions À AddedIdealUpPositions Statistical analysis was performed by generating 100 similar sized random gene lists, which were then evaluated with the same strategy to obtain the average value and the standard deviation. These values were then used to determine the significance level based on the calculated Z-score as described above. The values for 88 random genes are depicted in the Figures and other gene numbers led to very similar mean and standard deviations of the controls. Only for very small gene lists (below 20 genes) the standard deviation was getting larger, which impacts the significance of the UpRegScore for small clusters, such as Daf-21down_3 and Daf-21down_4.

GO enrichment analysis
GO Enrichment Analysis was performed with the help of the webserver at http:// geneontology.org/page/go-enrichment-analysis [39,40]. Here we uploaded gene lists from individual clusters and the enrichment analysis was performed in the category "Biological Process" without the Bonferroni Correction. Results were ranked according to their p-values and the highest ranking terms, which showed significant enrichment (marked with "+" in the enrichment column and p-values<0.001), were included in the bar chart Figures with their respective enrichment factor as y-axis value.

qRT-PCR for skr-5 and act-1
Nematodes were harvested 2.5 days after placement of L1 larvae on RNAi plates. Isolation of total RNA from the nematodes was performed with the SV Total RNA Isolation System (Promega Corporation) following the manufacturer's instructions. Primers used in the qRT-PCR experiments were designed with Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primerblast/). The final primer sequences were: act-1_fwd (AATCCAAGAGAGGTATCCTTA), act-1_rev (GATGGCGACATACATGGCT), skr-5_fwd (AATTGGTGCTGGCAGCCAC), skr-5_rev (GTTACC CAAGTTGAAAACGGCACG). qRT-PCR samples were prepared with the Brilliant III SYBR1 Master Mix Kit (Agilent). The PCR program on the Agilent Mx3000p was selected as described in the user manual. After the PCR reaction, the amplified DNA was tested to confirm the purity and the correct length of the product. This was done by on-chip gel electrophoresis with the Agilent 2100 Bioanalyzer system. We used the expression values only in cases, where we could obtain a specific product at the right size. The experiments were performed in triplicates and averaged.

Hsp90-downregulation leads to accumulation of chaperones and depletion of certain lectins
We first used mass spectrometry to investigate changes in the nematode's proteome upon Hsp90-depletion. To perform a mass spectrometric comparison between Hsp90-depleted and control nematodes, the daf-21 dsRNA-expressing bacterial strain was grown on plates containing solely 15 NH 4 SO 4 as nitrogen source. The control population was grown on NGM-plates containing 14 NH 4 SO 4 . RNAi against Hsp90 leads to slower development and sterility. Nematodes were harvested when they approached L4 larval stage with Hsp90-depleted nematodes being slightly smaller at the moment of harvesting (Fig 1A). After washing the nematodes off the respective plates the worms were lysed and equal amounts of the two isotope-specific samples were combined. The samples were then subjected to mass spectrometric analysis performing a direct comparison of the identified peptides in a single mass spectrometry experiment. To remove any impact the nitrogen source may have, the experiment was then repeated with exchanged nitrogen sources ( 15 NH 4 SO 4 for control, 14 NH 4 SO 4 for Hsp90-depletion) so that altogether two datasets were obtained to describe this system. About 400 proteins could be quantified and we compared their relative levels. Proteins upregulated after Hsp90-depletion in both experiments included many known heat-shock proteins and other proteins associated with the chaperone system ( Fig 1B, S1 Table). Similar to previous reports, these proteins were the strongest upregulated group, indicating that the induction of the heat-shock response is a significant event after depletion of Hsp90 [8]. As such, the small heat-shock proteins HSP-16.2, HSP-16.1 and HSP-16.48 are highly overexpressed (log 2 of about 4). STI-1, UNC-45 and other cochaperons of Hsp90 also showed increased protein levels, as does the protein folding machine CCT/TriC (log 2 of about 1, Fig  1B). Additional proteins found to be enriched in Hsp90-depleted nematodes were UNC-15, Y94H6A.10, Y41C4A.11 and F20D1.3. Interestingly, most of these differentially expressed proteins, including Y94H6A.10 and Y41C4A.11, can be clustered into a network using information from a transcriptional coexpression database ( Fig 1C, left side). This implies that these proteins are commonly found to be regulated in the same direction. Testing the upregulated proteins for enrichment of GO-Terms also shows a strong enrichment in "protein folding" and "unfolded protein binding" (Fig 1C, right side), as expected from the nature of the included proteins.
Several proteins also showed reduced levels in Hsp90-depleted nematodes ( Fig 1B, S2 Table). This group included several lectin-like genes (such as CLEC-4, CLEC-65, CLEC-66, CLEC-83) and aspartate proteases (such as ASP-2, ASP-5 and ASP-6). Successful knock-down was confirmed by Hsp90/DAF-21 showing up in this group as well (log 2 of -2.8). Here most proteins can be clustered based on their known coexpression ( Fig 1D, left side) and GO enrichment analysis indicates that these proteins are usually coregulated as part of the innate immune response (Fig 1D, right side) [38,[45][46][47][48]. This is true for the C-type lectins and aspartate proteases. Furthermore many of the unassigned proteins in this group, such as K12H4.7, T19D12.4, F55G11.4, C49C8.5, ZK6.11 and ZK228.4 also are associated with the innate immune response according to Wormbase [49]. Thus C. elegans obviously responds to Hsp90 RNAi with specific proteins originating from defined transcriptional clusters related to the heat-shock response and the innate immune response.

Gene expression changes upon depletion of Hsp90 highlight the impact on developmental processes
To better define the connection between these clusters, we investigated the full transcriptional program under this condition by employing Affymetrix DNA microarrays. Three experiments were performed under conditions similar to the proteome analysis by harvesting the nematodes shortly before they reached fertility. We averaged the expression differences for the biological replicates and obtained the genes that were most strongly up-or downregulated from the 26,959 sampled transcripts (Table 1). To calculate transcriptional networks, we used the top 250 genes in each direction and clustered these genes with the help of public coexpression databases as described in Materials and Methods [26]. To test the significance of the obtained networks, we compared their connectivity with 100 networks of random genes. Indeed, in both Hsp90-dependent networks, the connectivity parameters, such as the percentage of hits (yellow arrow) and the Hsp90-depleted nematodes are slightly smaller at that stage. B) Proteomic changes based on isotope-labelled mass spectrometry. The two independent experiments are plotted in the two dimensions. Blue spots indicate genes, for which the identification scores based on the number of detected peptides imply a reliable quantification. C) Network of upregulated proteins and corresponding GO enrichment analysis. The network was generated based on hit-to-hit relationships in a coexpression database as described in the Materials and Methods section. D) Network of downregulated proteins and corresponding GO enrichment analysis. The network was constructed with identical settings as in C. https://doi.org/10.1371/journal.pone.0186386.g001 Influences of Hsp90-downregulation in nematodes Table 1. Averaged expression differences of strongly affected genes after daf-21 RNAi. The genes depicted are the most strongly affected genes. The three experiments were averaged to obtain the average expression change and the standard deviation. These genes are part of the network structures of Fig  2. The weaker hits are truncated in this table, but included in the networks in Fig 2. The assignment to the clusters is based on the cluster description in Fig 2. "0" implies that the gene is outside of all marked clusters, "+" connections imply that the gene is in-between marked clusters and "-" implies that this particular gene was not found in the used coexpression database. In the column "strain" all currently available fluorescent reporter strains for these hits are shown, but strains in brackets were not used in our study. included (S1A Fig.) and the overall number of connections (S1B Fig) are much higher than for random genes (748 and 952 connections versus 12 and 16, p<0.001). This shows that the hits of this response appear commonly coexpressed. We further tested whether up-or downregulated genes outside the 250 hits can be successfully predicted from the coexpression database, a feature which shows the completeness of the transcriptional response. We quantified this ability with a score between -100 (worst possible prediction for hit genes 250-300) and 100 (perfect prediction for hit genes 250-300). The determined score is higher than 90 for both directions, while being close to 0 for the similar sized random gene groups (S1C Fig), showing that almost all predicted genes behaved as expected (p-value < 0.001). This confirms that both directions of our microarray experiment contain genes that are commonly coexpressed based on the used coexpression database [28]. This is not only true for the averaged hit lists but also for each of the biological replicates (S1C Fig), which confirms the high data quality of each individual experiment. The predictive ability evident from the high CoRegScores can only be explained if the Hsp90-depletion induced response shares many similarities with the hundreds of microarray experiments which were used to calculate the coexpression correlations in the modSEEK database [28]. We then visualized the connected hits as clustered networks (Fig 2A and 2B). The genes upregulated after Hsp90-RNAi arrange in two clusters (Daf-21up_1 and Daf-21up_2), which are minimally connected by coexpression relationships (Fig 2A). One of the clusters (Daf-21up_1) connects the genes of the heat-shock response (hsp-16.2 and hsp-70) previously shown to be induced after RNAi against Hsp90 [5,8,15] with 61 other genes. According to GO-term enrichment analysis and confirmed by individual searches, this cluster contains many genes related to the innate immune response (Fig 2A). It also includes the most strongly upregulated genes arf-1.1 (60-fold), Y41C4A.11 (50-fold), Y47H10A.5 (32-fold) and C50F7.5 (8-fold) ( Table 1). The large network of Daf-21up_2 instead originates from the slower development of the Hsp90-depleted nematodes. It contains genes involved in larval development and collagen and cuticle formation. Based on GO-term enrichment the genes from this cluster are mostly related to the molding cycle, cuticle development, and also to the development of individual tissues like pharynx and intestine (Fig 2A). We then tested whether these two clusters are upregulated in all three replicates alike. Indeed, in all replicates, the innate immune/ heat-shock response cluster Daf-21up_1 is strongly upregulated (UpRegScores of 86, 98 and 100). While the upregulation of Daf-21up_2 is slightly weaker in the first experiment (UpReg-Score of 70, 97 and 97, respectively), it is still clearly significant when compared to 100 random gene sets (UpRegScore: 3.7 ± 6.2, p<0.001) (Fig 3A). Thus the three biological replicates generate similar responses on the clustered level and this correlation is also evident from the same top hits being upregulated in all three experiments (S2- S4 Figs). Nevertheless, the biological replicates differ in the strength of their transcriptional responses, which may be caused by a variable RNAi-induced knock-down among replicates or slight differences of the developmental state at the harvesting moment.

Upregulated Gen Transgene Strain Average STD
Downregulated genes are involved in germline development and innate immune response.
We also visualized the clusters of genes downregulated after Hsp90-depletion (Fig 2B). Here two main clusters can be separated (Daf-21down_1 and Daf-21down_2) based on the database-retrieved coexpression relationships. Contrary to the upregulated dataset the connections between the two main clusters are much stronger. This implies that they may not be independently regulated. Two very small clusters (Daf-21down_3 and Daf-21down_4) are formed separately. Only the gene cyp-13A6 is not connected to any partner gene within the clustered network ( Fig 2B). Also for this side of the response, the three replicates show striking similarities and the UpRegScores for all clusters are between -54 and -99 compared to 3.7 ± 6 for random genes (p<0.001). Again the first replicate is showing a weaker response. According to GO-enrichment analysis and based on single-gene queries, the two main clusters contain genes that are involved in the development of different tissues from pharynx and epithelial cells to embryo development (Daf-21down_1 and Daf-21down_2). This also is evident from Wormbase information, where many genes in these clusters are associated with processes in the germ line (C39D10.7, Y62H9A.4, Y62H9A.6, ZC373.2 and C17G1.2), the spermatheca (clec-223, ZK813.3 and F55B11.3), or in the early embryo prior to gastrulation (pes-2.1, T04G9.7, vet-6 and F55B11.2). Thus the main part of the lower expressed transcriptional response reflects the slower developmental progress of the Hsp90 RNAi-treated nematodes and the stalled development of its gonads and embryos. In the first replicate where the emphasis of the response is shifted towards the immune response, highly inducible genes are primarily seen in the Daf-21down_4 cluster. Some of these genes, like the gene dod-24, have been described as DAF-16-regulated genes related to the innate immune response [50] and responsive to stress conditions [51]. Likewise, the dod-24 connected genes H20E11.3 and C32H11.4 are part of innate immune responses [47] along with the gene clec-209, which was found as top-ranking gene in innate immune responses against Vibrio cholerae [36]. Thus in the downregulated clusters, genes that are related to the innate immune response are set apart from the many genes that are connected due to their coexpression during the nematode's developmental process.
Expression changes after Hsp90-depletion overlap with described immune responses.
We next aimed at identifying conditions where similar expression patterns had been observed. We realized that several of the genes most prominent in our study had also been prominent in other microarray experiments, including a very recent study on the innate immune response [52]. We thus compared our response to these microarray experiments which were performed to study the effects of a depletion of SAMS-1 and SBP-1 [33,52]. We analyzed the publicly available data for RNAi against sams-1 and tested the similarity to Hsp90-depletion by plotting the sams-1-related expression values onto the networks obtained from our microarray data. To quantify the overlap in the respective clusters, we calculated the UpRegScore for each of our clusters ( Fig 3B, S5 Fig). Clearly the cluster Daf-21up_1 is regulated similarly in these two experiments (UpRegScore of 80, p<0.001) while the other clusters do not show an overlapping response (S5 Fig). The same is true for the response to depletion of sbp-1 (UpRegScore of 84, p< 0.001, Fig  3B, S6 Fig). Thus Hsp90-depletion induces significant parts of the innate immune response in a manner similar to other genes previously connected to the normal regulation of this response, like sams-1 and sbp-1 [33]. Furthermore, the exclusive induction of cluster Daf-21up_1 in the sams-1 and sbp-1 RNAi experiments confirms that Daf-21up_1 and Daf-21up_2 are indeed regulated entirely independent from each other. To directly test the overlap of Daf-21up_1 with the pathogen-induced immune response, we evaluated the response of our clusters during immune responses against Vibrio cholerae and Pseudomonas aeruginosa [36,37] (Fig 3B, S7 and S8 Figs). Indeed many of the genes of cluster Daf-21up_1 are also induced in response to these pathogenic bacteria (UpRegScore of 41, p<0.001 for V.cholerae and UpRegScore of 44, p<0.001 for P. aerginosa) but the overlap in this cluster is not as high as in the sams-1 and sbp-1 RNAi-induced responses. Again no significant similarity can be observed in the other clusters. Thus one part of the differentially expressed genes clearly shows that RNAi against Hsp90 strongly induces the innate immune response. The second set of differentially regulated genes instead shows the changes to the development of the nematode and also the impaired development of its embryos in an independently regulated response.

Hsp90-depletion induced changes to the immune response are evident in intestinal tissues
Having identified genes which are differentially affected after RNAi against Hsp90, we aimed at testing some of the main hits in fluorescent reporter strains. This also helps to define the tissues in which Hsp90-related changes can occur. Previous experiments had successfully visualized the induction of the heat-shock response in body wall muscle cells and intestinal cells upon Hsp90 knock-down [5,8]. For this aim, we initially obtained publicly available reporter strains for genes that are prominent in our immune response clusters. We aimed at investigating the strongest affected genes and preferred translational GFP-fusions over transcriptional reporters. The first strain contains a GFP-fused version of SKR-5. SKR-5 is a Skp1-homolog protein [53] expressed in intestinal tissues and upregulated 6-30 fold (log 2 : 3.8 ± 1.1) in our microarray experiments. We subjected these nematodes to RNAi against Hsp90 and control RNAi and then compared the fluorescence intensity and localization to see whether differential regulation can be observed in response to Hsp90-depletion. We see GFP-SKR-5-expression in some control treated SHK207 nematodes during later larval stages in posterior parts of the intestine as previously reported [44]. The Hsp90-depleted SHK207 nematodes instead show much stronger fluorescence in the entire intestinal tissue with most nematodes now being fluorescent in all cells of the intestinal tube (Fig 4A). The upregulation of this gene also is evident in qRT-PCR assays, where a three-fold upregulation can be observed (Fig 4B). We next tested a transcriptional clec-60 reporter strain. CLEC-60 is a C-type lectin known to be regulated as part of the innate immune response. It is expressed in larval intestine and in adults in the Int8 and Int9 cells [54,55]. In control RNAi experiments, we observed the reported behavior. Upon RNAi against Hsp90, the fluorescence is strongly increased (Fig 4C). To exclude differences that may arise from the slow growth of Hsp90-RNAi treated nematodes, we tested if earlier larval stages of the control-treated nematodes show an expression with similar intensity to that of Hsp90-dsRNA treated nematodes. However, this was not the case for these reporter strains (data not shown), implying that the observed upregulation in intestinal tissues reflects the consequences of the Hsp90-depletion. We then tested a transcriptional reporter strain containing dod-24p::GFP. DOD-24 is an epoxide hydrolase connected to the innate immune response [51]. This strain shows high expression throughout the whole intestine under normal growth conditions (Fig 4D, left side). After treatment with Hsp90-RNAi, the expression in the Int3 to Int7 intestinal cells is strongly reduced while it is retained to a weaker extent in the very anterior and posterior intestinal cells ( Fig  4D, right side). Thus, in these three reporter strains which change their expression as part of the nematode's immune response, the intestinal expression of fluorescent markers correlates well with the expression analysis, confirming the influence of Hsp90-RNAi on this immune response pathway.
Genes from developmental clusters are affected in the nematode's gonads and embryos.
We then tested strains for hits that cluster in the early developmental networks. Here we initially picked a translational reporter strain for SEPA-1, a protein, which mediates autophagy of P-granules in early embryos [41,56]. The expression of this gene product, which is degraded once embryos reach 100 cells, is decreased in our microarray data upon Hsp90 knock-down. Specifically SEPA-1 (Fig 5A) can be observed in the syncytial gonad arms and in the embryos prior to and just after passage through the spermathecum. We tested the influence of Hsp90depletion on the localization of this protein and its expression. Indeed, we see reduced fluorescence of this stage-specific protein and incomplete development of gonad arms containing this protein. Moreover, the embryos did not show SEPA-1 containing structures after passage through the spermathecum, implying that fertilized oocytes are not generated correctly. We then tested a translational reporter for the caveolin homolog CAV-1 [42] from the Daf-21down_2 cluster. This protein is required during the first cell division and then quickly disappears from the embryo. Oocytes approaching the spermatheca in Hsp90-depleted worms contain CAV-1, but are reduced in size and no fertilized oocytes expressing CAV-1 can be seen after the passage through the spermatheca (Fig 5B). However, in few cases, very brightly fluorescent oocytes are found containing high amounts of CAV-1 without evident structural organization (data not shown). These likely represent the previously observed endomitotic oocytes [4]. These changes show that Hsp90 depletion significantly impacts the developing oocyte prior to fertilization. We finally tested the expression of the bZIP-like transcription factor ZIP-8 [43], which is also present in the Daf-21down_2 cluster. Based on the fluorescent reporter strain BC12422, pzip-8 is activated in early embryos and is also present at later stages in the adult intestine. We did not observe pzip-8::GFP expression in oocytes of Hsp90-depleted nematodes, while it was clearly visible in Influences of Hsp90-downregulation in nematodes L4440-control treated nematodes (Fig 5C). Thus the reduced expression levels within this cluster point to the failure of Hsp90-depleted nematodes to develop functional oocytes.
We finally tested a reporter strain of the downregulated gene pud-1.1. While pud-1.1 not included in the main clusters of the downregulated network, it is loosely connected to cluster Daf-21down_1. PUD-1.1 is one of two identical PUD-1 proteins of unknown function that are known to be differentially regulated during aging and development [57]. In control nematodes, Influences of Hsp90-downregulation in nematodes GFP-PUD-1.1 is found in the intestine at later larval stages [58]. We see a decrease of fluorescence in the intestine of Hsp90 RNAi-treated nematodes, implying that the upregulation in developing adults does not occur in Hsp90-depleted nematodes (Fig 5D).
These strains collectively confirm that processes during the development of oocytes are affected by depletion of Hsp90. This apparently happens in addition to the changes in intestinal tissues, where heat-shock response and innate immune response are observed.

Discussion
Depletion of Hsp90 by RNA interference leads to various morphological and transcriptional changes in C. elegans. The combined phenotype of the Hsp90 knock-down includes developmental changes to gonad and vulval structures, induction of the heat-shock response, changes to muscle ultrastructure, and, as observed here, the induction of the innate immune response. Based on the function of Hsp90 as a regulator of signaling kinases and transcription factors, these diverse changes could be caused by multiple Hsp90-dependent processes. Here we used whole genome expression data to define pathways that are linked to Hsp90-depletion.

Separate transcriptional pathways are induced by RNAi against Hsp90
To understand the induction of separate pathways, it is important to organize the differentially expressed genes into their respective expression clusters and then attempt guesses on transcription factors for the coexpressed gene groups. We recently implemented a gene clustering approach that was based on the use of genome-wide coexpression data. This yielded information-rich networks in yeast that clearly separated independent clusters and thus divided hit lists into separate transcriptional units [26,27]. Potential transcription factors could then be derived from lists of target genes as available from YEASTRACT [26,59]. Here we apply the same clustering strategy to the multicellular nematode C. elegans and like in yeast information-rich networks are obtained with all tested parameters proving a highly significant network and cluster formation. In this study, the C. elegans response to Hsp90 knock-down can be separated into multiple distinct clusters that can be further defined as independent based on the very low interconnection-numbers. This is supported by independent microarray experiments, which trigger only part of the Hsp90-depletion response. GO-term enrichment and individual searches further show that this clustering approach also separates the genes according to distinct functional processes. The reporter strains used in this study validate our organization strategy since genes from the same cluster behave similar regarding their tissue-specific expression and their fluorescence changes after Hsp90 RNAi. By this approach the response to Hsp90-depletion can be separated into one part that reflects a strongly induced innate immune response and into another part that reflects the slower and incomplete development of the nematodes and their reproductive structures.

HSPs are upregulated on the proteome and transcriptome level
The induction of the heat shock response after Hsp90-depletion has been observed in reporter strains before [5,8] and our study confirms this on the protein and mRNA levels. The transcriptional networks overlap with the proteomic response in particular for the heat-shock proteins HSP-16.1, HSP-16.2 and HSP-70. On the protein level the wider chaperone network is elevated, including many of the Hsp90 cofactors such as CDC-37, UNC-45, STI-1, ZC395.10 and FKB-6. Even though hsp-16.2 and hsp-70 are included in the transcriptional networks, the expression changes for most chaperones are barely evident on the mRNA level. This also is true for Hsp90 itself, whose mRNA is only weakly suppressed at the moment of harvesting (log 2 (DiffExp) = -0.15). This may imply that the chaperone system in the harvested nematodes has already built up a new balance at the beginning of phenotype-development. The heatshock response certainly is not the strongest response on the transcriptional level at the moment of harvesting as other transcriptional clusters are more prominent among the top 250 genes.

Hsp90-depletion induces the innate immune response in intestinal cells
The induction of the immune response is much more evident at the harvesting stage. Here, the upregulation of some genes, such as Y41C4A.11 (log 2 of 5.69) and Y94H6A.10 (log 2 of 2.42) from cluster daf-21Up_1 (innate immune response/heat-shock response) can also be observed on the proteomic level (log 2 of 3.06 and log 2 of 3.31, respectively). While Y41C4A.11 is a distant C. elegans homolog of the lipopolysaccharide-induced TNF-factor alpha, the function of Y94H6A.10 is unknown but its general coexpression with the heat-shock response and the innate immune response is evident in all experiments of this study. In microarrays and mass spectrometry approaches applied in this study, immune response genes are found to be both upregulated and downregulated. As such, several C-type lectins of the innate immune response are reduced after RNAi-treatment against Hsp90 [38,[45][46][47][48]. In general, most of the proteins found at reduced levels in our mass spectrometry experiments share a functional or coexpression connection to components of the innate immune response. This result holds true at the transcriptomic level, as many upregulated genes from cluster daf-21Up_1 (C08E8.4, Y47H10A.5, B0024.4, F15B9.6) can also be found upregulated in other cases where the innate immune response is activated. This is similar to the immune response after contact with Vibrio cholerae [36] or Pseudomonas aeruginosa [37] and for the response to RNAi against sams-1 and sbp-1 [52], two genes that are apparently able to suppress the innate immune response like daf-21.
As expected, the expression in reporter strains for skr-5, clec-60 and dod-24 is most affected by Hsp90 RNAi in the intestine, the place where most reactions of the immune response take place. It is interesting that several prominent hits of this response group are regulated by the transcription factor DAF-16, a known regulator of the immune response [60][61][62]. skr-5 and dod-24 are direct target genes of DAF-16, along with hsp-16.2, fbxa-163, C08E8.4, F13D11.4, and Y105C5A.12 [44,51], which also are among our top induced genes. Given that DAF-16 target genes are well described, we could test whether certain clusters of our networks are strongly enriched for DAF-16 targets, which would suggest that DAF-16 activity is influenced in the Hsp90 depleted nematodes. To this end, we employed a recently published genomewide ranking of DAF-16 targets, which assigned 1663 Class I targets (upregulated by  and 1733 Class II targets (downregulated by DAF-16 and its co-regulator PQM-1) [44,63]. Using this resource, we determined the ranking number for each clustered gene (Fig 6A, S9  Fig) and tested whether any of our clusters show significant DAF-16 correlation. It is evident that cluster daf-21Up_1 in particular contains genes which are regulated by DAF-16 (Fig 6B  left plot, S9 Fig). All told, 71% of the genes in this cluster are assigned as DAF-16 targets, which is significantly enriched over the 18% expected by normal distribution. Clusters daf-21Down_1 and daf-21Down_4 also show enrichment of DAF-16 targets, but they are formed from genes repressed by DAF-16/PQM-1 (Fig 6C, S9 Fig). These data suggest that the DAF-16 pathway is highly active after Hsp90 RNAi. It is interesting in this respect that the induced DAF-16/PQM-1 response differs from the published target ranking, as a group of 12 genes typically repressed by DAF-16/PQM-1 are now activated in daf-21Up_1. This may suggest that the balanced regulation between DAF-16 and its transcription suppressing co-regulator PQM-1 could be altered by Hsp90 RNAi, leading to the induction of the innate immune response as observed here. Germline-specific markers confirm effects of Hsp90-depletion on events before fertilization Germline development is disrupted in Hsp90-depleted nematodes. This is evident from the sporadic failure of Hsp90-depleted nematodes to reach fertility and develop correct gonadal or vulval structures [4]. The germline-specific changes were not visible on the proteomic level, likely because most of these proteins are not among the 400 proteins we could quantify. However the transcriptional response does show very consistent genome-wide changes. In cluster Daf-21down_1 and Daf-21down_2, we find many genes which have a relationship to germline and embryo development. This matches the disruptive influence of Hsp90 RNAi on development. It nevertheless has to be kept in mind that GFP-expression from genomic integrated plasmid arrays may not entirely correlate with the native context of this protein. Thus, a more detailed analysis of the affected cell types would need to involve detection methods for endogenous mRNA or protein, like fluorescence in situ hybridization or antibody staining. Despite this the tendency is obvious that most tested target genes from this cluster are coexpressed within the gonads and developing embryos. Thus, reporters for the promoter activity and protein localization of cav-1, zip-8 and sepa-1 will be helpful to get information on the delayed processes during development. How these changes originate is still enigmatic but the identification of defined developmental alterations related to the observed phenotype could serve as a starting point to study the failures during development after Hsp90-knock down in more detail.
Could the fertility problems also be responsible for the immune response? Several studies have been published about the influence of germline development on the immune response and they suggest that disruption to early germline development can induce the immune response genes by a soma-protective pathway [60,64,65]. It remains unclear whether the incomplete embryo development as observed after Hsp90 RNAi can induce the soma-protective immune response. Only in sporadic cases where the gonad arms are obviously shortened by the Hsp90 knock-down are the meiotic germ cells also affected. In our networks, the genes overexpressed after Hsp90 RNAi share very few coexpression relationships with the large gene cluster related to the growth defects. This would imply that the upregulation of these immune response genes usually is not chronologically coupled to the growth defects but is independently regulated instead. To gain more information on this, we individually observed genes which should be upregulated if the soma-protective pathway is active, which would connect the reduction in mitotic germ cells to the activation of the DAF-16 and ATF-7 transcription factors. An important gene, upregulated as part of this pathway is irg-7/F40F4.6, which influences the DAF-16 regulated longevity genes and the PMK-1/ATF-7 regulated immunity genes [66]. We see very limited upregulation of irg-7 in our microarray data (log 2 (DiffExp) = 0.51) and we do not see a systematic influence on PMK-1/ATF-7 targets [38, 67-69] (F08G5.6: log 2 (-DiffExp) = 0.14, C09H5.2: log 2 (DiffExp) = 1.03, C17H12.8: log 2 (DiffExp) = 0.25, F56D6.2: log 2 (DiffExp) = -0.14, T24B8.5: log 2 (DiffExp) = -0.99, K08D8.5: log 2 (DiffExp) = -0.08). Given that these genes should be activated via irg-7 this may imply that the soma-protective route via irg-7 is not triggered extensively after Hsp90-depletion. While these relationships need to be investigated in further detail, the presented variations provide a starting point to disentangle the complex relationships following the depletion of Hsp90 in C. elegans. each cluster the position of the genes is marked by a point at the respective ranking position. The top 1663 genes represent targets considered upregulated by DAF-16; the bottom 1733 genes represent targets downregulated by DAF-16 according to this ranking. B) Histograms for the two clusters daf-21Up_1 and daf-21Up_2, where the bar on the left side roughly represents the 1663 Class I targets and the bar on the right side the 1733 Class II targets. C) Histograms for the clusters downregulated after Hsp90 RNAi. https://doi.org/10.1371/journal.pone.0186386.g006 Influences of Hsp90-downregulation in nematodes Supporting information S1 Fig. Microarray analysis of Hsp90-depleted nematodes. Parameters obtained from the cluster analysis for the top 250 differentially expressed hits are shown for each of the three RNAi experiments, which evaluate the significance of the networks (first replicate: column 1 and 2, second replicate: column 3 and 4, third replicate: column 5 and 6) and the average (column 7 and 8). A) Percentage of hits included in the network for all three experiments. As comparison the results for random gene lists (columns 9 and 10) are depicted. A hit gene is counted as included in the network if a single coexpression connection to any other hit gene is found in the coexpression database. Most of the connections in the random networks are only isolated hit-to-hit pairs. B) Average counts of connections in the network per hit. Two genes will yield many connection counts, if this pair is marked as coexpressed several times in the database. C) CoRegScore for each network to quantify the predictive strength of the network. The score is calculated as described in the Materials and Methods section.  [63]. Red nodes indicate Class I targets, which are upregulated by DAF-16. The shadings imply the ranking within this group of 1663 genes, with the most intense red tone marking the strongest DAF-16 targets. Blue indicates Class II targets in this ranking. These genes are downregulated by DAF-16 together with the transcription factor PQM-1. The most intense blue tone indicates the strongest class II target. Genes, which could not be found in this ranking are omitted, the 15890 genes not considered DAF- 16