Caenorhabditis elegans hub genes that respond to amyloid beta are homologs of genes involved in human Alzheimer’s disease

The prominent characteristic of Alzheimer’s disease (AD) is the accumulation of amyloid beta (Abeta) proteins in the form of plaques that cause molecular and cellular alterations in the brain. Due to the paucity of brain samples of early-stage Abeta aggregation, animal models have been developed to study early events in AD. Caenorhabditis elegans is a genetically tractable animal model for AD. Here, we used transcriptomic data, network-based protein-protein interactions and weighted gene co-expression network analysis (WGCNA), to detect modules and their gene ontology in response to Abeta aggregation in C. elegans. Additionally, hub genes and their orthologues in human and mouse were identified to study their relation to AD. We also found several transcription factors (TFs) responding to Abeta accumulation. Our results show that Abeta expression in C. elegans relates to general processes such as molting cycle, locomotion, and larval development plus AD-associated processes, including protein phosphorylation, and G-protein coupled receptor-regulated pathways. We reveal that many hub genes and TFs including ttbk-2, daf-16, and unc-49 have human and mouse orthologues that are directly or potentially associated with AD and neural development. In conclusion, using systems biology we identified important genes and biological processes in C. elegans that respond to Abeta aggregation, which could be used as potential diagnostic or therapeutic targets. In addition, because of evolutionary relationship to AD in human, we suggest that C. elegans is a useful model for studying early molecular events in AD.


Introduction
Extracellular senile plaques and neurofibrillary tangles of tau protein are the main hallmarks of Alzheimer's disease (AD). Senile plaques are composed of amyloid beta (Abeta) peptides with neurotoxic features that lead to synaptic dysfunction, connectome disruption, and neural death. In familial AD, Abeta deposition is related to mutations in either Abeta precursor protein (APP) or catalytic subunits of γ-secretase, presenilin-1 (PS1) or presenilin-2 (PS2) [1]. γ-PLOS ONE | https://doi.org/10.1371/journal.pone.0219486 July 10, 2019 1 / 28 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 secretase and β-secretase are membrane-bound proteases that cleave APP to Abeta. γ-secretase cleaves the C-terminal fragment of APP and generates multiple types of Abeta, based on the position of cleavage [2]. Two dominant forms of Abeta40 and Abeta42 constitute approximately 80-90% and 5-10% of total Abeta in normal cells, respectively. However, Abeta42 is the major molecule in the formation of Abeta plaques, which is observed in one-third of plaques [3]. Abeta peptides are involved in multiple physiological processes, including regulation of synaptic activity, transcription, neuronal survival, processing of APP, and antioxidant activity [4][5][6]. Removal of Abeta from cells is achieved via neprilysin (insulin-degrading enzyme) followed by lysosomal degradation (REF). It can also be transported into blood vessels and carried away for degradation [1]. Despite the importance of these peptides in AD, it is not possible to monitor the accumulation trends and its impact on human neural cells. Therefore, the role of Abeta and its accumulation is analyzed in animal model systems of AD. Transgenic modifications have generated models of AD in animals including primates [7], mouse [8], rat [7], fruit fly [9], and nematode (Caenorhabditis elegans) [10]. Although C. elegans is a convenient model for neurodegenerative disease, it differs from humans in lacking Abeta and β-secretase [10]. These animals do not exhibit human pathological and behavioral symptoms of AD, but they can be used for studying molecular mechanisms related to Abeta aggregation. Using transgenic techniques, several C. elegans strains have been generated to express human Abeta, where Abeta precipitates in muscle and leads to paralysis [10]. To study early responses to Abeta aggregation and to identify proteins and pathways affected by Abeta, model systems have been introduced that harbor temperature-dependent induction of Abeta expression [11]. In such models, the gene expression profiles can be studied using highthroughput techniques at the onset of Abeta accumulation as described in the upcoming sections.
Transcriptome data can also be analyzed using network approaches to determine modules of genes or proteins, molecular mechanisms, interactions between members of a system, and the most important proteins in a pathway that could be used as potential candidates for further experiments [12]. Weighted Gene Co-expression Network Analysis (WGCNA) is a robust method, widely used to detect important members (genes) and modules of genes that permit downstream functional analysis [13]. WGCNA estimates the pair-wise correlation between genes, using expression values, and identifies groups of genes with a similar pattern of expression [14]. Another informative network analysis method is the Protein-Protein Interaction (PPI) network approach. Proteins can function in complexes, therefore studying these interactions is effective in identifying key members in networks. Information to construct such networks originates from either experimental findings or in silico predictions [15].
Here, we have applied a systems biology approach, using data from transgenic C. elegans models of AD that synthesize Abeta protein, to find early molecular responses and potentially important genes and proteins in the response to Abeta. Several studies have used AD biomarkers that show the onset and progress of the disease, however, they are still not able to determine molecular modifications in brain tissue at early stages [16]. We have detected orthologous genes in worm, mouse and human suggesting that similar pathways are affected by Abeta expression in these evolutionarily distant species.

Microarray data description
In order to study the early effects of Abeta on C. elegans at the systems level, we searched the Gene Expression Omnibus (GEO) database of National Center for Biotechnology Research (www.ncbi.nlm.nih.gov/gds) using "amyloid beta" keyword and filtered the outcome by selecting C. elegans as the organism. We found that among all resulting datasets only one study under accession number GSE65851, originally published by Hassan and colleagues 2015, would be useful for pursuing our goal (having a sufficient sample number and experimental design within early stages and including a successive time-course). In this study, Abeta toxicity was investigated in C. elegans, using a strain (CL4176) that expresses human Abeta42 in body wall muscle. Toxicity was examined at 4-hour intervals from time-point 0 to time-point 6 (20 hours) [11]. This dataset contains 47 microarrays, 15 of which are related to Abeta toxicity and the remaining are related to GFP expression, as the negative control, or GFP-degron to study the impact of Abeta protein aggregation on cells [11]. Briefly, they performed microarray analysis, gene ontology RNA interference, and paralysis assays on samples collected from several time-points with 4-hour intervals. The sample collection covered time points between T0 and T20, where severe paralysis begins. Through these time points the process of early response to Abeta could be studied. Here, we performed network-based analysis and gene ontology to find relationships between genes and potential hub genes (Fig 1). Furthermore, to independently validate our findings, two additional microarray datasets from the human brain were also analyzed. These two datasets were GSE12685 that encompass 6 AD and 8 control samples from the original study by Williams and colleagues 2009 [17], and GSE28146 containing 7 incipient AD and 8 control arrays published by Blalock et al. 2011 [18].

Data preparation and gene expression analysis
To obtain DEGs (Differentially expressed genes) between different time-points and experimental conditions, GEO2R tool (www.ncbi.nlm.nih.gov/geo/geo2r/) was used, which is embedded in NCBI. For the Abeta study in worms, two sets of analysis were performed. First, the transcriptome profiles of Abeta-manipulated samples were compared to that of GFP expressing samples. While in the second analysis, the expressions of genes were studied in samples over consecutive time-points. Besides, at each time point, the gene expression of GFP samples were removed from that of Abeta accumulation to exclude developmentally linked DEGs. Significantly affected genes were identified by applying two filters; first, selecting only those genes with a |log2 FC| of 0.6 between two conditions and then a p-value � 0.005 (for the time-course study) or p-value � 0.01 (for the comparing made between Abeta and GFP samples). Similarly, in datasets that were used for independent validation, genes harboring a pvalue � 0.05 and showing ±0.6 changes in Log2FC were considered as differentially expressed genes (DEGs). Next, all of the probes without annotation were removed. For duplicate probes, those with the lowest p-values were included in the analysis. To prepare the required highquality input matrix for the WGCNA package, one needs to determine the distribution of the samples and data quality for network analysis. To achieve this, we first performed principal component analysis (PCA) and filtered outliers. To further eliminate bias and noise from the data, the matrix of gene expression normalized by log2 transformation and probes with average log2 expression value less than 0.6 (through all arrays) were discarded. Based on the coefficient of variations, only 0.9 percentile of the genes were selected as input for WGCNA.

Gene coexpression analysis via WGCNA
WGCNA analysis was performed according to Langfelder and Horvath, 2008 [14]. Applying R package WGCNA, we constructed an "unsigned" network for all arrays using "softpower = 12", and "deep-split = 2". The remaining default thresholds were used with no further changes. Module eigengenes were obtained and their relationships with different time-points were determined. At each time-point, highly related modules were selected to determine the correlation between module membership and gene significance. Then modules with the highest correlation were identified and used for further analysis. Networks were visualized and analyzed using either Cytoscape 3.4.0 [19] or Gephi 0.9.2 [20].

PPIs network among differentially expressed gene products
Several "undirected" PPI networks were constructed using DEGs and interactions between proteins were extracted from the STRING database [21]. This was achieved by submitting DEG lists to the database and selecting information for C. elegans. Interactions were filtered for being "experimentally validated", "information obtained from other databases", and "coexpression", and meet a medium confidence threshold of 0.4. PPI networks were analyzed and visualized by Cytoscape 3.4.0 [19] and Gephi 0.9.2 [20]. To detect protein complexes as subnetworks, we used ClusterONE plugin in Cytoscape software [22], considering the following settings: minimum size: 5; minimum density: Auto; Edge weight: string confidence value. All the presented modules had p-values � 0.005. The top 10 percent of nodes with the highest weighted degree were selected as hub genes in the PPI network.

Transcription factor analysis
The list of TFs in C. elegans was retrieved from Reece-Hoyes and colleagues study [23]. They identified a set of 934 TFs by computational analysis and manual curation. Here, using hierarchical clustering and correlation method, we clustered the TFs in our datasets based on their log2 normalized expression values through all arrays. Additionally, differentially expressed TFs (DE-TFs) in at least two stages or overexpressing DE-TFs (see the results section) were selected and clustered. The clustering was performed by Heatmap.3 in R.

Gene ontology and enrichment map construction
We used the Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 for gene ontology analysis [24]. For each gene list, we have selected up to 5 biological process terms with the lowest p-value (� 0.05). For constructing an enrichment map of biological processes we have used Enrichment map 3.0.0 plugin of Cytoscape, using default settings [25].

Homology analysis and literature review
To identify human and mouse orthologous genes for the hub genes in C. elegans, "tblastn" was used on the NCBI server (blast.ncbi.nlm.nih.gov/Blast.cgi). First, FASTA format of the longest protein sequences were fetched for the selected genes and submitted to "tblastn" and blasted against Homo sapiens and Mus musculus using default parameters. Then, the top hits with Evalues below 1e-10 were obtained and considered as the human or mouse counterparts of the input genes.
A literature review was performed by searching "Google" (www.google.com), "Google Scholar" (www.scholar.google.com), and "PubMed" (www.ncbi.nlm.nih.gov/pubmed) with keywords, including "[the gene name]" with either "Alzheimer", "Abeta", or "Amyloid Beta" terms. This approach helped us to find potential publications directly related to the genes of interest. Next, the resources were manually checked to determine the role of the given gene in AD (if any). In doing so, we have determined possible changes in gene expression, methylation, and protein level alterations in addition to single nucleotide polymorphism in the gene of interest in AD. Consequently, a relationship between the gene of interest and AD or any other disorders in neuron activity was established (direct or indirect). We also conducted similar approaches for gene aliases.

Data selection
Here, we have analyzed transcriptome data produced by Hassan and colleagues in which they used a transgenic C. elegans model expressing Abeta in body wall muscle [11]. It is known that C. elegans is a suitable model to study neurodegenerative diseases including AD [26].

PCA and DEGs analysis
We performed PCA analysis to detect the distribution pattern and to identify any potential outliers of the microarray samples containing gene expression data of Abeta-expressing worms at different time-points [12]. The samples were collected at 4-hour intervals from timepoint 0 to 6, thereafter denoted as T0, T4, T8, T12, T16, and T20. Our analysis show that the samples of each time-point are clustered with minor overlaps between T8 and T12 (Fig 2A). Notably, only one repeat in T0, T4, and T20 failed to group with their counterpart samples. As the samples were grouped as expected (no outlier was found), all samples were included in the remainder of our analysis. Next, DEGs were identified using GEO2R, a tool embedded in NCBI, that uses R package Limma. As this tool requires at least two samples in each treatment, we were not able to compare the gene expression pattern of time-points T0 and T4. Analyzing the other time-points, we detected the highest number of DEGs in the comparison between time-points T16 and T12 (with 2156 DEGs), which is in agreement with the PCA result. While the lowest number of DEGs (273 DEGs) was detected for the comparison made for the timepoint T8 vs. T4. The highest differences between up-and down-regulated DEGs observed in the comparison of T8 vs. T4 (with 241 and 32 genes being up-and down-regulated, respectively), and in the results obtained for the T20 vs. T16 comparison (where 13 and 398 genes were up-and down-regulated, respectively). In two other comparisons, namely T16 vs. T12 and T12 vs. T8, an almost equal number of up-and down-regulated genes were observed ( Fig  2B). Venn diagram representation was used to show the common DEGs among different comparisons ( Fig 2C). The highest number of common DEGs were found between the comparisons made for time-points T16 vs. T12 stage and T12 vs. T8, with 184 DEGs. Only 13 DEGs were found in all comparisons including T16 vs. T12, T12 vs. T8, and T8 vs. T4. Interestingly, no genes were present in all time-points, indicating time-point specific expression trends.
Gene ontology enrichment analysis showed that different biological processes are associated with the different stages of transition ( Fig 2D). Metabolic processes (6.3% of DEGs involved in this processes) and cell shape regulation (1.9% of DEGs) are the most affected processes at the transition from T4 to T8 (named as stage one). While, during the transition from T8 to T12 (stage two), processes including nematode larval development (15.3%), locomotion (12.7%), and body morphogenesis (6.7%) are represented by the highest number of DEGs. Protein phosphorylation was dominant process at the transition from T12 to T16 (stage three), however, cell shape regulation (2.3% of DEGs) was also observed. Annotation of DEGs in the T16 to T20 (stage four) transition indicated that G-protein coupled receptor signaling pathway (7.9% of DEGs), and detection of chemical stimuli involved in sensing (3.9% of DEGs) were affected the most. Furthermore, enrichment map analysis was performed to detect connections between these processes and to identify the most central processes at each transition ( Fig 2E). Maximum deregulation of genes occurred at transition from T12 to T16. c) A Venn diagram for the comparison of the DEGs among time-points after Abeta expression. d) Affected biological processes in each time-point comparison. The numbers close to the bars indicate the percentage of DEGs that are involved in each process. e) Enrichment map of biological processes, in which nodes are biological processes and edges show common genes between two nodes. Larger nodes represent those processes that contain more DEGs. Main process in each network is indicated by red star. There was no enrichment network for the T8 vs. T4 comparison. In PCA graph large dots are average of all the sample of that group.
In the second stage, locomotion and nematode larval development were the most important processes. While, protein phosphorylation and G-protein coupled receptor signaling pathways were the most central processes in the third and fourth stages, respectively. No enrichment map was detected for the first stage.
When Abeta and GFP samples were compared (Fig 3), we found a lower number of DEGs compared to the time-points comparisons. Here, the highest number of affected genes were at time-point T16, where 825 DEGs were detected (Fig 3A and 3B). Gene ontology indicates that proteolysis, metabolic processes, and neuropeptide signaling pathways are important biological processes associated with these DEGs (Fig 3C).

WGCNA construction and module detection
We have detected modules and highly connected genes related to each time-point by applying WGCNA to the list of gene expression values. Hierarchical clustering grouped the time-points together, with some exceptions. No outliers were detected in the samples, which is in agreement with the PCA results. Soft-threshold analysis revealed that soft-power 12 would be the best value for module detection (Fig 5A). Following network construction and module detection, module-time-point relationships were also investigated (Fig 5B and 5C). In total, 8 out of 21 detected modules were mostly related to three time-points T8, T16, and particularly T12. To show exact relationships with these time-points, module membership vs. gene significance correlation analysis was performed for the most related modules (Fig 5D). We selected the top 8 modules that were highly correlated with at least with one of the time-points. It should be mentioned that some modules, in spite of being highly correlated, did not have enough members to pass the significance threshold for gene ontology analysis. Altogether, modules named  Table of top three biological processes in the modules from as "antique white 4", "brown 4", "cyan", "dark grey", "dark olive green 2", "dark violet", "purple", and "violet" were identified and gene ontology enrichment was performed for them. Several important biological processes were affected by the expression of Abeta protein at networks. Note-the biological processes with P-value less than 0.05 were selected. In the network, modules are indicated by numbers. Nodes and labels in larger size indicate higher centrality-based degree factor. Thicker edges indicate higher weight. There was no network for time-points T8 and T12.
https://doi.org/10.1371/journal.pone.0219486.g003  Table contains information regarding the top three biological processes found for the genes in the modules detected in each network. The biological processes with pvalue less than 0.05 were selected. In the network, modules are indicated by numbers. Red and blue colors of nodes represent up-and down-regulation, respectively. Nodes and labels in larger size indicate higher centrality based on degree of connectivity factor. Thicker edges indicate higher weight. https://doi.org/10.1371/journal.pone.0219486.g004 Similarities between responses to Abeta accumulation in worm and human different time-points, including G-protein coupled receptor signaling pathway, transport, regulation of transcription, meiotic nuclear division, embryo development, locomotion, and reproduction (Table 1). In order to reduce the number of the modules and increase the significance of the results, the modules with similar gene ontology and lower correlation with the time-points were removed and only those with highest significant in GO were considered for further analysis. Finally, six modules including "antique white 4", "brown 4", "cyan", "dark grey", "dark olive green 2", and "purple" were selected for network analysis and to detect the Similarities between responses to Abeta accumulation in worm and human hub genes in response to Abeta protein. Interestingly, we found that many of the hub genes in the modules identified as "antiquewhite4", "cyan", and "purple" were differentially expressed. Expression patterns of these hub DEGs were similar among all modules, where the genes were down-regulated in modules labeled as "antiquewhite4", and "purple", whereas they were upregulated in the module named "cyan" (Fig 6). Similarities between responses to Abeta accumulation in worm and human Identified modules were labeled as a) "antiquewhite4", b) "brown4", c) "cyan", d) "darkgrey", e) "darkolivegreen2", and f) "purple".

Modulation of transcription factor expression in response to abeta accumulation
By comparing the 582 TFs proposed by [23] against the expression matrix prepared in the current study, 540 TFs were detected to be present in the transcriptome of C. elegans in response to Abeta accumulation at different stages. Cluster analysis based on similarity method has identified several clusters of TFs across the different time-points. However, only two of these clusters showed an altered pattern of expression according to the progression through timepoints ( Fig 7A). Many of the TFs were not differentially expressed, therefore we selected only the DE-TFs that are present at least in two comparisons. Six DE-TFs, including uaf-2, daf-16, lin-29, egl-27, unc-62, and hlh-10 were detected and labeled as DE-TFs cluster. Interestingly the expression patterns for these DE-TFs followed similar trends in the comparisons, except for the hlh-10 that showed slightly different expression trends to that of other DE-TFs ( Fig 7B).
In addition, to find common TFs that express similarly in response to accumulation of Abeta protein, a cluster was identified and labeled as over-expressing TFs, which contains only those DE-TFs that show accumulative trends through consecutive stages. These up-regulated DE-TFs were selected and re-clustered to identify those TFs that respond to the accumulation of Abeta protein. Over-expression was observed for the TFs, including tbx-34, pos-1, oma-1, mex-6, ceh-60, mex-5, fkh-6, cey-2, nhr-234, lin-11, and cey-3 through successive time-points, therefore these TFs could be considered as the major regulators of gene regulation upon increases in the Abeta expression ( Fig 7C).

Specific responses to Abeta accumulation during development
The time point comparisons we conducted in the previous sections do not rule out the possible impact of developmental stages on the observed perturbation of gene expression. To identify specific factors which solely respond to Abeta accumulation (and not developmental stages), we conducted similar analysis but instead of comparing each time point with the previous one, at each time-point we compared gene expression of Abeta-GFP expressing samples with those of GFP expressing counterpart (S1 and S2 Figs). Through WGCNA analysis we identified 3 modules specifically responding to Abeta accumulation and not GFP expression (Fig 8). The GO analysis of the members of these modules including "Sienna3", "greenyellow" and "yellowgreen" are provided in S3 Fig Biological processes such as detection of chemicals, locomotion and metabolic process are among the most affected pathways. Interestingly, the specific responses to Abeta accumulation were evident at 12 and 16 hours after induction. We also extracted and analyzed the subnetworks of the involved modules (S4 Fig). The most important TFs and genes in the comparisons were compiled and presented in Figs 9 and 10, which also include literature review for the role of each entry.

Homology analysis, gene relationships to AD and validation of results
Human and mouse orthologues for the hub genes and DE-TFs were identified through Blast analysis. Numerous worm genes had at least one orthologue in the human or mouse genome (Figs 9 and 10). Regarding WGCNA, only differentially expressed hub genes were selected and their orthologues were identified. Results showed that 65% (26 out of 40) of the hub genes detected by WGCNA in worm had an orthologue in either mouse or human or both Red and blue colors of nodes represent up-and down-regulation, respectively. Grey nodes represent no DEGs for that gene. Nodes and labels in larger size indicate higher centrality in the network based on degree factor. Thicker edges indicate higher weight in terms of network analysis. Only fifty percent of edges are shown. https://doi.org/10.1371/journal.pone.0219486.g006 organisms. Most orthologous genes were detected in modules labeled as "antiquewhite4" and "purple". In contrast module "brown4" contained the lowest number of orthologues (Fig 9). Hub genes detected in the PPI networks of C. elegans had the lowest percentage of identified orthologues, where we could identify orthologues for only 11 out of 32 genes. However, when only DE-TFs were blasted against the human and mouse genomes, almost all of these TFs that affected by Abeta accumulation in worm had at least one orthologue in the other organisms, except for hlh-10 (Fig 10).
In the next step, we reviewed the literature to assess if the orthologous genes, detected in human and mouse, are linked to AD. Additionally, expression patterns of orthologues were also investigated in two other independent datasets available for incipient AD. We found that 58% of the orthologues (15 out of 26) for the genes detected as a hub in the WGCNA analysis had a direct link to AD. In addition, we established indirect links to AD for the remainder of the orthologous genes in that analysis. Furthermore, 26% of the orthologous genes were directly linked to Abeta perturbation (Fig 9). Regarding the orthologous genes detected for the hub genes in the C. elegans PPI networks, 7 out of 11 had a direct link to AD, while, 6 out of 11 of these orthologous genes had a potential linkage to AD. Regarding the transcription factors, orthologues of DE-TFs were also associated with AD. We found that while gene expression perturbations for 56% of identified orthologues were reported in AD, 3 of these DE-TFs were Similarities between responses to Abeta accumulation in worm and human directly linked to Abeta accumulation. Table 2 provides a convenient view of the hub genes, with information related to the function, cellular compartment, and the cell type or organ that the genes are expressed in, in addition to the function of their human orthologous genes. Interestingly, some of the genes are expressed in neural cells, however, they are also expressed in other tissues.
We validated our analysis with datasets related to incipient AD. After removing outlier arrays, as totally separated samples from the clusters (GSM318211, and GSM697312 from GSE12685 and GSE28146, respectively), DEGs, common DEGs, and biological process were detected (Fig 11). Analyzing these datasets showed variation in the transcriptome profile between different studies in the incipient stage of AD. Processes related to transport and synaptic transmission were common in both datasets. Comparing the list of orthologous genes with the list of genes obtained from these two datasets identified that several of human orthologues of the Abeta-responding genes in the worm undergo gene expression modifications in these datasets, particularly GSE12685. However, when we analyzed the trends of gene expression, different trends of expression were observed (Figs 9 and 10). Most orthologous genes  Similarities between responses to Abeta accumulation in worm and human  obtained from WGCNA and DE-TFs analysis of worm data were found to be differentially expressed in the human AD datasets (Fig 10).

Discussion
Owing to the high prevalence and mortality rate, AD has been the center of attention for many neuroscience researchers. Although numerous genes, proteins, and molecular processes related to the pathogenicity of AD have been identified, the key genes and proteins involved in the onset of the disease have remained mostly elusive [1,27]. The main obstacle in this type of research is inaccessibility of brain tissue samples at the early stages of AD, which is the asymptomatic period of the disease. As an alternative approach, early molecular alterations in AD can be examined in animal models that develop disease symptoms. Among these animal models, C. elegans expressing human Abeta protein enables detection of changes in the expression of genes in response to Abeta aggregation. In spite of the genome similarities between human and worm, there is a large evolutionary and structurally distance between these animals that brings into question the reliability of associations between the genes obtained from worm to that of human disease.
In order to address this issue, we used two network approaches, WGCNA and PPI, to detect important genes and pathways that respond to the accumulation of Abeta in C. elegans. PPI networks were constructed using DEGs from analysis at different temporal stages of Abeta-GFP expression compared to GFP-expressing control samples. Gene modules and hub nodes related to each stage were obtained from WGCNA, which uses original expression values rather than differentially expressed values. Generally, we found several biological processes relevant to AD in both analyses of WGCNA or PPI. Transcriptome alterations in response to Abeta accumulation were most evident at the later time-points (T12 and T16). However, the alteration of gene expression decreases from T12 onwards, suggesting that Abeta protein accumulation might cause gradual alterations to gene expression and exerts its impacts rather early in the progression of AD. This observation could explain why the onset of AD could happen many years before recognizing symptoms [1].
We could link many of the identified biological processes in the C. elegans response to Abeta protein directly to AD in human, including metabolic processes [28], protein phosphorylation [29], G-protein coupled receptors [30], meiotic chromosome segregation and cell cycle [31], and regulation of transcription [32]. These similarities reflect possible evolutionary conservation in the molecular mechanisms of responses to Abeta in both worms and humans. There were also alterations in other processes including the molting cycle, collagen, and cuticulin-based cuticle, olfactory behavior, body morphogenesis, locomotion, hermaphrodite genitalia development, and nematode larval development. Additionally, similar to what the original study reported regarding detection of heat shock proteins in response to Abeta, we detected genes responsible for heat shock responses at stage 2 (T12 vs. T08) (Figs 2 and 3) [11]. Interestingly, similar to what we found here, a previous study indicated that overexpression of APP-related protein (APL-1) in neurons disrupts olfactory behaviors in C. elegans [33].
C. elegans is different from humans in many anatomical, molecular and physiological characteristics [34]. Most notably, C. elegans lacks β-secretase or Abeta expressing genes in its status in two datasets from human incipient Alzheimer's disease, the direct relationship with Alzheimer's disease, any potential role in Alzheimer's or association with Abeta accumulation were also included. For each gene, higher intensity in color shows higher expression. NF indicates homologous genes could not be found. � TFs responding to Abeta accumulation through development, �� TFs uniquely responding to Abeta accumulation when compared to GFP sample at the same developmental stage. Brown and green bars indicate trends of expression in response to GFP and Abeta accumulation, respectively. Darker color indicates more expression. https://doi.org/10.1371/journal.pone.0219486.g010 Similarities between responses to Abeta accumulation in worm and human Table 2. Encoded protein, gene ontology and anatomical location of hub genes in worm and biological processes of their orthologues in human. All the Information are extracted from Uniprot (https://uniprot.org) and Wormbase (https://wormbase.org) databases. Anatomical location indicates cell/organ/regions that the genes are expressed in. Dash lines means no information has been found.

In worm
In human  genome [10]. Additionally, in these two organisms, Abeta aggregation occurs in different tissues; in muscles of C. elegans model (in this study CL4176) compared to the brain in human [11]. In spite of these differences, our systems-level observations indicate that Abeta aggregation affect somewhat similar processes in both systems. Responses to Abeta aggregation are expected to be exerted via hub genes in networks constructed for these systems. Three modules ("antiquewhite4", "cyan", and "purple"), identified by WGCNA, contained the most affected, differentially expressed hub genes (Fig 6). We observed synchronized trends of expression among these hub genes as well as hub proteins that were identified in the PPI network analysis. Next, the relationship between these hub genes and proteins to AD or Abeta was investigated by mining previously published studies. Although there are publications describing molecular mechanisms of C. elegans models for AD [35,36], as we expected, there are few publications studying AD or Abeta accumulation in C. elegans. We only found two hub genes in C. elegans that have already been previously reported to be associated to Abeta accumulation, including mec-12 that encodes alpha tubulin subunit [37], and daf-16 that mediates proteostasis to tolerate β-amyloid toxicity [38]. Therefore, for the first time, we have established some degree of association between these hub genes and proteins with Abeta accumulation in C. elegans.
Following Abeta accumulation, the expression level of multiple TFs may be affected, and as a result gene regulatory networks in cells would be perturbed. We found that several transcription factors are deregulated in response to Abeta accumulation (Fig 7). This finding is in agreement with a previous study on the role of daf-16 in response to Abeta aggregation in C. elegans [38]. Additionally, daf-16 and its orthologous FOXO genes are major regulators of aging and stress responses [39,40]. It has been shown that mutation in daf-16 increases dopaminergic neurodegeneration, however, DAF-16 is not the only neuroprotective agent [41]. In our study, daf-16 was not continuously over-expressed and the highest expression occurred at T12 which could be in response to severe muscle paralysis. The importance of two other TFs, EGL-27 and NUC-62, in postponing aging and increase in C. elegans lifespan have also been reported [42,43]. The important roles of these TFs in the aging process and their alteration by Abeta aggregation, suggest potentially crucial roles for them and their human orthologous TFs in AD. Interestingly, the literature review for the human orthologous of the over-expressing TFs showed a high percentage of association with AD (Fig 10). RERE is orthologue of EGL-27 that is involved in AD has been detected by bioinformatical methods and gene expression analysis [44]. ZNF362 and PBX1 are human orthologue of LIN-29 and CEH-60 that are also deregulated in post-mortem AD patient brain samples [45,46]. Deletion of Nr2f6 (Ear2), the mouse orthologous of NHR-234, leads to defects in early memory and learning in the mouse model for AD, APP/PS1 [47]. This point suggests that these over-expressing TFs are important in responding to Abeta aggregation and start their activity later than other TFs.
By searching for orthologues in human or mouse genome we found that many of the hub genes have orthologous counterparts in these organisms that are directly or potentially involved in AD (References are included in Figs 9 and 10). Generally, in WGCNA, the hub genes in modules related to transport and determination of lifespan ("antiquwhite4") and nematode larval development "purple" had the highest number of orthologs in the human and mouse genome. Interestingly, gene ontology of the orthologous genes revealed that they are present in pathways including neuroactive ligand-receptor interaction plus glutamatergic and dopaminergic synapses. Among these genes, we found some orthologues with a remarkable role in AD. For example, TTBK2 (the C. elegans C39H7.1 orthologue) over-expression reduces phosphorylation of tau protein in AD [48]. While CCKAR (ckr-2 in C. elegans) has been identified as a biomarker of AD [49]. Another important gene, HNRNPA1 the human orthologue for sqd-1, is involved in alternative splicing of APP genes that could be targeted for reducing senile plaque formation [50]. Polymorphism in the PPP3CA gene (orthologue of tax-6 in in C. elegans) has been reported in AD patients, indicating a role for this gene in the disease [51]. APP endocytosis and APP alpha-secretase cleavage is under SNX33 regulation, which is an endocytic protein orthologue of lst-4 in worms [52]. Treusch and colleagues studied toxicity of Abeta in yeast and found several genes that modify Abeta toxicity. Interestingly, they found orthologues of these genes are modifiers of Abeta toxicity in glutamatergic neurons of C. elegans and primary rat cortical neurons [53]. Two of our hub genes, unc-26 and lst-4, are present in this list. Remarkably, their orthologues in mammalians are also linked to AD [54,55]. Out of the genes discussed above with a potential role in AD, C39H7.1, lst-4, and daf-16 appear to specifically respond directly to the accumulation of Abeta (Figs 9 and 10), while others could be affected either by Abeta or developmental progression. However, we could not rule out a combinational impact of both these factors on gene expression changes. On the other hand, pos-1, nhr-234, mex-5 and mex-6 are all overexpressed in response to Abeta accumulation and all have been directly or indirectly linked to AD [47,[56][57][58][59]. This shows the efficiency of our approach to detect the most important factors in response to Abeta accumulation.
It has been shown that gene regulatory networks are conserved across metazoa [60]. As expected, altered TFs in C. elegans had orthologues in human. Interestingly, many of the human orthologous TFs are directly associated with AD or other neurodegenerative disorders. This suggests that conserved regulatory networks are involved in responses to Abeta aggregation and AD pathogenesis. We reasoned that deregulation of the human orthologous genes might be altered during early AD. To examine this, we have selected two datasets from incipient AD in human. We found deregulation in a number of these genes, however, opposite trends of expression were observed in different datasets. Unexpectedly, we found that two datasets are also different in DEGs and expression patterns (Fig 11). Many of these genes are previously reported, indicating the involvement of these homologous genes in responses to Abeta accumulation in both human and worm (PMID of the references are included in Figs 9 and 10).
In conclusion, using systems biology we have identified important genes and biological processes in C. elegans that respond to Abeta aggregation. These results could be useful for examining the molecular mechanisms involved in Abeta aggregation as a main cause of AD, which may be applied as further diagnostic or therapeutic targets. Additionally, we suggest C. elegans could be considered as a useful model for studying early molecular events in AD because of the evolutionary relationships to AD in humans documented in this study.