hsp-90 and unc-45 depletion induce characteristic transcriptional signatures in coexpression cliques of C. elegans

Nematode development is characterized by progression through several larval stages. Thousands of genes were found in large scale RNAi-experiments to block this development at certain steps, two of which target the molecular chaperone HSP-90 and its cofactor UNC-45. Aiming to define the cause of arrest, we here investigate the status of nematodes after treatment with RNAi against hsp-90 and unc-45 by employing an in-depth transcriptional analysis of the arrested larvae. To identify misregulated transcriptional units, we calculate and validate genome-wide coexpression cliques covering the entire nematode genome. We define 307 coexpression cliques and more than half of these can be related to organismal functions by GO-term enrichment, phenotype enrichment or tissue enrichment analysis. Importantly, hsp-90 and unc-45 RNAi induce or repress many of these cliques in a coordinated manner, and then several specifically regulated cliques are observed. To map the developmental state of the arrested nematodes we define the expression behaviour of each of the cliques during development from embryo to adult nematode. hsp-90 RNAi can be seen to arrest development close to the L4 larval stage with further deviations in daf-16 regulated genes. unc-45 RNAi instead leads to arrested development at young adult stage prior to the programmatic downregulation of sperm-cell specific genes. In both cases processes can be defined to be misregulated upon depletion of the respective chaperone. With most of the defined gene cliques showing concerted behaviour at some stage of development from embryo to late adult, the “clique map” together with the clique-specific GO-terms, tissue and phenotype assignments will be a valuable tool in understanding concerted responses on the genome-wide level in Caenorhabditis elegans.


hsp-90 and unc-45 depletion induce characteristic transcriptional signatures in coexpression cliques of C. elegans
Lukas Schmauder & Klaus Richter * Nematode development is characterized by progression through several larval stages. Thousands of genes were found in large scale RNAi-experiments to block this development at certain steps, two of which target the molecular chaperone HSP-90 and its cofactor UNC-45. Aiming to define the cause of arrest, we here investigate the status of nematodes after treatment with RNAi against hsp-90 and unc-45 by employing an in-depth transcriptional analysis of the arrested larvae. To identify misregulated transcriptional units, we calculate and validate genome-wide coexpression cliques covering the entire nematode genome. We define 307 coexpression cliques and more than half of these can be related to organismal functions by GO-term enrichment, phenotype enrichment or tissue enrichment analysis. Importantly, hsp-90 and unc-45 RNAi induce or repress many of these cliques in a coordinated manner, and then several specifically regulated cliques are observed. To map the developmental state of the arrested nematodes we define the expression behaviour of each of the cliques during development from embryo to adult nematode. hsp-90 RNAi can be seen to arrest development close to the L4 larval stage with further deviations in daf-16 regulated genes. unc-45 RNAi instead leads to arrested development at young adult stage prior to the programmatic downregulation of sperm-cell specific genes. In both cases processes can be defined to be misregulated upon depletion of the respective chaperone. With most of the defined gene cliques showing concerted behaviour at some stage of development from embryo to late adult, the "clique map" together with the clique-specific GO-terms, tissue and phenotype assignments will be a valuable tool in understanding concerted responses on the genome-wide level in Caenorhabditis elegans.
Nematode-development is a highly complex process that is defined by temporal and spatial events in different tissues and cell types. Therefore simultaneous events are occurring in this process with chronological timing to enable the highly reproducible development program.
HSP-90 (DAF-21) is a molecular chaperone, crucial for the development of vulva, gonads and oocyte maturation as well as ensuring longevity of C. elegans [1][2][3] . It is an indispensable protein, activating and regulating many clients, for example protein kinases, and transcription factors, such as steroid receptors [4][5][6] . Inhibition of HSP-90, by either RNAi or specific compounds, therefore has the potential to interfere with several pathways. RNAi arrests the nematode development and reduces motility in later larval stages 7,8 . Prominent responses induced after hsp-90/daf-21 RNAi include the heat-shock response, which is known to be suppressed by HSP-90 in most organisms 1,9,10 . Other affected responses are potentially regulated in a more organism-specific manner, like the innate immune response, which is coupled to the heat-shock response in nematodes 11,12 . Interestingly, both of these responses are also dependent on the developmental state of the nematode, with the heat-shock response being barely inducible in early larvae and also in adult aging nematodes 11 . The reason for these correlations is unclear, but it could be supported by assigning genes clearly to individual responses, so that the common principles and regulatory patterns become obvious.
The HSP-90 cofactor UNC-45 participates in the muscle-specific functions of HSP-90. Invertebrates possess a single unc-45 gene, which is expressed in muscle cells, where UNC-45 performs HSP-90-dependet folding of the myosin motor domain. It further is expressed in non-muscle tissues of early embryos [13][14][15]  www.nature.com/scientificreports/ or blue tone was defined in Cytoscape (https:// cytos cape. org/) 23 . In cases where responses were weak, like both unc-45 experiments, the scale was adjusted to reach from log2 < − 0.25 to log2 > + 0. 25. This analysis leads to information on most cliques as to whether they are induced or repressed with statistical significance as described before 17 . This method to evaluate nematode array data will be implemented for public use in the clusterEX.de webserver, which currently has this functionality only for yeast arrays. Correlation analysis between different samples was made by plotting the cliques' expression values against each other and obtaining the coefficient of determination R 2 for the regression line. If R 2 was closer to 1, the correlation between the two sets was considered to be better. These results were compared to correlations on the gene level in cases where identical array types were utilized.
Analysis of microarray data on development. To define moments of clique relevance during development, time points from developmental series were used to determine a transcriptional status for each clique in the map. In many cases, cliques react to developmental steps as concerted units resulting in a non-random distribution of up-and downregulated hits throughout the 307 cliques. To cover several larval states, three published GPL200 series were obtained from the GEO repository. These represent a time course of early development with data from embryo, L1 and L4 larvae (GSE6547 24 ) and a time course describing the aging process with time points at L4 larvae, and adults at day6 and day15 of development (GSE21784 25 ). Lastly, a time-course was included describing different stages of larval development, composed of L3, L3-lethargus, L4, L4-lethargus and young adult (GSE46291 26 ). Expression values were obtained from the normalized data table containing all public GPL200 experiments (see above).

Results
A genome-wide coexpression clique map for the nematode C. elegans. To obtain transcriptional units influenced by hsp-90 and unc-45 RNAi-treated nematodes, we first generated gene cliques that are coregulated in C. elegans, in which each of the 22,620 genes is assigned to exactly one clique. We had used the same procedure before to generate a coexpression clique map for S. cerevisiae 17 . Based on the same stepwise procedure, we grouped every gene reported on standard microarrays of the GPL200 platform into one coexpression clique of at least five genes. The procedure resulted in 307 coexpression cliques, which were visualized in Cytoscape to generate the "coexpression clique map" for C. elegans. We set out to test, whether these 307 coexpression cliques are gene groups with a high level of functional similarity, as it was observed for the yeast clique map before 17 . Therefore, we investigated all cliques by GO-term enrichment analysis. 220 of the 307 cliques show a GO-term enrichment with a p-value lower than 1e-3, 172 showed less than 1e-4, and 148 of the 307 cliques showed p-values of less than 1e-5 (Best results in Table 1). This is far better than 307 random cliques, which yielded these p-values 18 times, 2 times and zero times. We also found significant enrichments employing phenotype enrichment analysis (PEA 22 ) and tissue enrichment analysis (TEA 21 ) with 145, 106 and 81 cliques being enriched for the same phenotype (20, 3 and 0 times in random cliques) and 45, 37 or 29 cliques being enriched for tissuespecific expression in the three p-value categories (0, 0, and 0 times in random cliques). The values also are far better than cliques composed of random genes. In this way, roughly two thirds of the 307 coexpression cliques were assigned either a function, a related phenotype or a tissue-specific expression with acceptable significance criteria of below 1e-4 (Table 1).

hsp-90 RNAi affects embryo development and induces stress responses.
Having confirmed that the "clique map" of coexpressed genes also holds information on functional, phenotypic and tissue-specific signatures, we set out to investigating the transcriptional response of hsp-90 RNAi-treated nematodes. We previously had analysed these microarray data based on the Top250 differential regulated genes obtained from three experiments 11 . hsp-90 depleted nematodes showed sterility, incomplete development of gonad arms and the formation of endomitotic oocytes 7,27 . Development is mostly blocked at a later larval stage. TAC analysis revealed many genes with substantially different expression levels and showed the strongest response in the experiment 1 (P152), while the experiment 2 (A966) and 3 (P062) showed a weaker response (Fig. 1a). Gene expression changes had implied the induction of the heat-shock response and the innate immune response in analyses before 11 , but due to the focus on only 500 of the 22,620 genes measured with this array type, information from the many weaker affected genes could not be considered in this study 11 .
The genome-wide gene expression cliques as defined here, instead allow visualization and analysis of all values. Significance analysis showed enrichment of up-and downregulation in many cliques for experiment 1 (Fig. 1b), but also for the other experiments 2 and 3 (Supplemental Fig. 1b and c). Indeed almost half of the cliques respond to the RNAi-treatment with a concerted response of their genes (Best cliques in Table 2). We first determined the upregulated cliques: these contain the clique col-138-col-49, which is holding genes related to the "structural constituent of cuticle" and the clique abu-7-abu-8_22491 related to the "response to unfolded protein". The largest upregulated clique, containing 209 genes, is mlt-9_22518-F33D4.6_14044 ("cuticle development" with phenotype of "molt variant" and localized in the "embryo hypodermis") and other cliques related to cuticle formation, including R12E2.14_75-R12E2.15, col-117-col-167_1015, col-146-col-133, col-128-cdh-10_9234 (enriched in "peptidase activity") and sqt-2-dpy-9. Cliques which are strongly upregulated also in experiment 2, include the cliques related to the "immune system response" C10C5.2-Y58A7A.3 and K08D10.9-F46A8.1. Based on the assignment of GO-terms, phenotypes and tissues, the largest and strongest downregulated cliques ( Table 2) represent "embryo development" (T22D1.5-inx-14, enriched phenotype of "aneuploidy" and localized in "embryonic germline precursors"), "embryo development" (T24D1.3-egg-1, enriched phenotype of "polar body defective early embryo") and "reproduction" (puf-3-oma-2_18268, enriched phenotype "meiotic chromosome segregation variant" localized in the "germline_precursors") among many other cliques hinting at the stalled www.nature.com/scientificreports/ gonad development in agreement with the sterility phenotype observed. Comparing the three experiments a weak correlation can be found between experiment 2 and 1 (R 2 = 0.148) on a gene-by-gene level, which is increased, if cliques are compared (R 2 = 0.307, Fig. 1c). The same trend can be seen between experiments 3 and 1 correlating with R 2 = 0.622 for a gene-to-gene comparison, which increases to R 2 = 0.764, if cliques are compared (Fig. 1d).
Moreover the significance analysis employing 20 random cliques shows that the most strongly up-and downregulated cliques also are usually fulfilling the 1e-5 significance criterium in the compared experiments ( Fig. 1c and d, colored in red). www.nature.com/scientificreports/ unc-45 RNAi leads to delayed conclusion of sperm and vulva development. We then investigated the RNAi treatment against the HSP-90 cofactor unc-45 with the same approach. unc-45 RNAi-treatment leads to developmental disruptions and incomplete fertility at a more adult stage. To see, whether differences in the cliques can be observed we performed two independent RNAi-experiments with subsequent transcriptome analysis on DNA microarrays. Analysis with TAC showed a weaker response compared to the hsp-90 RNAi in both experiments (Fig. 2a, Supplemental Fig. 2a). This also was evident in the analysis of the 307 expression cliques, where the color scheme had to be adjusted to visualize the concerted reactions ( Fig. 2b and Supplemental Fig. 2b). Gene-gene comparisons showed a coefficient of determination of 0.15 between the experiments. When cliques were compared a coefficient of determination of 0.48 was obtained (Fig. 2c), confirming that also very weak responses can yield higher levels of repeatability by comparing matched groups of genes and not individual genes. Like with hsp-90 RNAi, specific cliques were found in all experiments to be significantly altered in their expression behavior. Upregulated are a few smaller cliques, like col-117-col-167_1015, msp-63-msp-33 and abu-7-abu-8_22491. These represent decisions to produce cuticle collagens, linker cell movement and induction of a response to heat. This also is reflected by the induction of hsp-16.2-F44E5.4_19238, a small clique containing the heat-shock proteins. This induction of the heat-shock response has been observed for unc-45 RNAi before 7 . Downregulated cliques represent the large cliques ZK1053.4-C08F1.6 (embryo development), ZC373.2-Y62H9A.6_1596 (cell membrane biogenesis) and sdz-10-fbxb-62 (L1 larval development). The weak differences, while being statistically significant for the described cliques, correlate with the mostly adult state of the nematode after unc-45 RNAi treatment and the observation that most developmental steps were performed, but the correct embryo development and the development of the vulva structure were affected 7 .
A hint towards the lacking germ line development may be derived from the misregulation of the cliques msp-36_msp-63 (linker-cell migration variant), which contains several genes related to sperm development and the downregulation of nspc-1_614-nspc-10_22525 (spermatogenesis variant). Thus, while vulva development and sperm development are stalled, certain features of the regulatory pathways are not deviating from the control nematodes and only later stages of the programmatic decision process show deviations that could explain the mismanaged development in the absence of unc-45.
Expression in developmental stages is altered in similarity to the RNAi-induced arrest. Having observed cliques with altered expression, we aimed at understanding, whether these expression changes are specific for one developmental transition occurring at the time point of arrest. We thus generated a time series of development ranging from embryo to late adult and compared the expression of all 307 cliques and in particular of those found relevant for unc-45 RNAi.
Striking differences were observed, when comparing the stages of each series (Supplemental Fig. 3), while differences between experiments of the same stage were small (Supplemental Figs. 4 and 5). Interestingly, also in these comparisons most of the isolated expression cliques showed coordinated expression differences, and also strong responses could be observed for the later developmental stages (Supplemental Fig. 6). In total more then 80% of the cliques show a statistically significant expression change during the development from embryo to 16 day adult and this also relates to most cliques found affected after unc-45 RNAi (Fig. 2d, 2 cliques and their development). While only few cliques were affected upon unc-45 RNAi treatment, hsp-90 RNAi is expected to yield a much stronger response.
Indeed a drop is observed in the expression of most upregulated cliques between L4 and day6 adult. In these cases the developmental delay may be the reason of the observed higher expression. A opposite pattern is observed for the downregulated gene cliques, with the exception of two cliques, which are not appropriately regulated: T05E12.6_12439-T05E12.6_12396 and gpd-3_977-aldo-1_21168, both of which appear to regulate metabolism.
Expression in developmental stages is altered in similarity to the hsp-90 RNAi-induced arrest. We next tested, whether also for the hsp-90 RNAi-treated nematodes developmental stages can be defined. The complexity of the differential expression between hsp-90 RNAi arrested nematodes and young adults allow to compare the obtained expression patterns with know patterns from larval development. We thus were interested to see, whether the full extent of the transcriptional changes can be explained by the observed developmental delay. Therefore we utilized publicly available microarray experiments on nematode development to help identify transcriptional units in the clique map that report on comparable steps during   www.nature.com/scientificreports/ development. We employed microarray data from three experimental series (Table 4) and initially compared developmental transitions, showing similarity to the differences we observe in the RNAi-treated nematodes. These comparisons were L3/young adult, L4/young adult and L4let/young adult (Fig. 3a-c) as investigated in GSE46288/GSE46289 28 . Clearly similarities can be observed between the hsp-90 RNAi treated nematodes and the L4 larvae, when each of them is compared to the young adult control. In fact, most of the cliques correlate in color and correlation analysis shows a coefficient of determination with these data of 0.4046, 0.5913 and 0.5915 (Fig. 3d). Based on these values, hsp-90 RNAi-arrested nematodes best correspond to a L4-larval like state. Only few clear differences can be observed compared to L4 or L4-lethargus, while several cliques deviate from L3-like state. Judged from the few differences to L4 state, it might be that the chronological timing of the events during development is misaligned in hsp-90 RNAi-arrested nematodes. We further investigated, whether the expression behavior matches the known expression behavior during development and aging. To this end the cliques found relevant for hsp-90 RNAi were investigated throughout development. While the strongly responsive upregulated cliques col-138_col-49, abu-7_abu-8, R12E2.14_75-R12E2.15, col-117_col-167, pqn-54_abu-9, col-146_col-133 already were found in context with the unc-45 RNAi-arrested state (see Fig. 2d), further cliques are identificed as significantly upregulated in hsp-90 RNAi treated nematodes (Fig. 3e). As for unc-45 RNAi before the downregulated cliques vit-2_vit-4, C46C2.5_W03F11, ZC373.2-Y62H9A.6_1596, ZK1053.4-C08F1.6 were identified amid fbxc-28_sdz-28 (see Fig. 2d), plus additional ones not found significant before (Fig. 3e). Interestingly, in few cases, the directionality of the RNAi-induced response does not correlate with the expected behaviour during the arrested L4 state. This is evident for the cliques C17H1.6-C17H1.13, C32F10.4-D1086.2 and C10C5.2-Y58A7A.3.

daf-16 target genes deviate from the developmental program in hsp-90
RNAi. We finally aimed at understanding the few cliques that deviate from the developmental progression and the explanation based on developmental delays. To this end we used the information gained previously that a fraction of the misregulated genes are daf-16 targets 11 . We tested, which of the cliques from the clique map contain daf-16 targets and then tested, whether those are regulated in coordance with developmental progress. Indeed, targets upregulated and suppressed by DAF-16 are enriched in several cliques (the 15 most prominent shown in Table 5, more information in Fig. 4a and b). Comparing the clusters identified in Eckl et al. (2017), with the current cliques we also observe a clear enrichment among several of the 307 cliques ( Table 5). As spectulated in Eckl et al., among the cluster "Up1" there are many genes, which are regulated by DAF-16, while cluster "Up2" does not enrich daf-16 targets (Table 5). Mapping all cliques onto the network developed in Eckl et al. the enrichment of these cliques in certain parts of the network becomes evident. For the downregulated genes, also DAF-16 enriching cliques are among those containing these genes 11 . Therefore, especially among the upregulated genes, cliques are present, which contain an elevated level of daf-16 target genes.
Interestingly, these cliques are upregulated despite their developmental program, which aims for downregulation. Thus, the presence of these cliques suggests a simultaneous modification to the dauer-program outside the developmental program after hsp-90 RNAi induced growth arrest.

Discussion
In this study, we analysed microarray data from C. elegans based on preformed coregulated expression cliques. This approach has been applied successfully in the yeast model organism by us 17 , but the applicability of this method to multicellular organisms has not been clear. We thus used the algorithms developed for yeast to also generate high quality coexpression cliques from nematode expression data and then validated them by GOterm enrichment, phenotype enrichment and tissue enrichment and by selective clique responses in individual microarray experiments. Based on our data from the developmental process of C. elegans, we believe that this analysis method could have broader use in the analysis of gene expression data from nematodes. This is evident from the correlated responses of cliques during nematode developmental transitions.
Recently also a different approach was reported to utilize genome-wide co-expression cliques for C. elegans 29 . In contrast to our appraoch, in this study an individual gene could be assigned to multiple cliques (on average 3) and also negative correlation was included. This makes the construction of a static clique map as used by us more difficult, but may include details missed by our approach. Both approaches will have their advantages. In the method described by us, we focus on the strongest connection and blank out those that might be secondary based on numbers, but still achieve very high levels of correlation with GO-terms, phenotypes and tissue specificity for most of the cliques.
One way to use the cliques could be by employing the popular GSEA platform 30 , where our cliques can be either used as a single input file covering the whole genome or as part of the global collection of gene sets. Another way to use the cliques can also be via the clusterEX.de webserver that we have set up and will further develop for the purpose of gene expression analysis based on known co-expression relationships. It therefore will be interesting to see, how further useful applications will be developed based on these predefined gene sets.

Integrating unc-45 into the developmental time line exposes distinct cliques for develeopmental stop.
We first analysed unc-45 depleted nematodes. In these nematodes, the depletion of unc-45 leads to developmental arrest and paralysis in almost adult animals. Here the comparison with the young adult nematode shows that certain cliques are misregulated and some of those cliques also represent developmental marker cliques as suggested by our evaluation procedure. These marker genes help to map the developmental status of the unc-45 depleted organisms. Clearly unc-45 depleted nematodes are close to N2 nematodes in this approach, but defined changes in certain genes help to map the events that did not unfold during development.
To evaluate the disruption of vulva development, we individually tested the genes transcriptionally regulated during this process and their specific regulation (Table 3) Regarding the germline, asb-2 is reduced (− 0.21, tars-1-AFFX-r2-3026-5_at) and the nspd-proteins are still upregulated together with msp-proteins (Fig. 2d 34 ), implying that sperm development is not completed yet, while the expression of the upstream regulators spn-4 and neg-1 35 is at the same level as in the normally developed adult. Also the regulators of msp-expression set-17 and csr-1 are expressed at the level of the control nematodes 36 , implying that sperm-development is almost finished 37,38 . Integrating hsp-90 into the time line data exposes defined clusters for developmental stop. We used this clique map to also analyse the depletion of hsp-90. While depletion of hsp-90 leads to developmental arrest and reduced motility in late larval stages, it also leads to defined transcriptional changes. To analyse the causes, we performed microarray experiments under wildtype conditions and under conditions, where the chaperone is depleted. Based on the clique analysis, it is obvious that certain developmental milestones are not reached yet in the HSP-90 depleted animals. Based on this analysis these nemtodes arrest in a late larval stage with additional misregulation of DAF-16 target genes. a.  www.nature.com/scientificreports/ Previously it had been observed that the Top300 genes from the hsp-90 RNAi analysis showed partial overlap with daf-16 regulated genes. We thus employed the gene-list from this previous study to identify the cliques, which now represent these genes. Indeed the correlation is fairly clear, with the cliques C17H1.6-C17H1.13, C32F10.4-D1086.2 and C10C5.2-Y58A7A.3 being mostly overlapping with the previous cluster1_up and the cliques col-138-col-49, R12E2.14_75-R12E2.15, mlt-9_22518-F33D4.6_14044 being mostly overlapping with the cluster2_up. Utilizing the ranked list of daf-16 target genes, we also determined which cliques most strongly are enriched in the Top750 and Bottom750 of this ranked list. These cliques are found mostly in cluster1_up confirming that the identification of this correlation also is visible from the clique map. Interestingly these cliques represent those that are differently regulated compared to the L4 larval stage. Thus the HSP-90 depletion leads to higher exression levles in a daf-16 regulated cluster (cluster1_up) and a daf-16 independent cluster (cluster2_up). With the daf-16 independent cluster containing mostly cliques related to larval development, apparently the depletion of HSP-90 induces both of these processes. Whether they are connected via secondary effects is unclear to date, especially as the developmental timing of DAF-16 activity is a well described phenomenon.

Scientific Reports
Thus, based on several clearly regulated marker cliques, hsp-90 arrested nematodes, like unc-45 arrested nematodes, can be positioned in respect to a developmental time axis.  Table 3     www.nature.com/scientificreports/ www.nature.com/scientificreports/