Co-expression analysis identifies neuro-inflammation as a driver of sensory neuron aging in Aplysia californica

Aging of the nervous system is typified by depressed metabolism, compromised proteostasis, and increased inflammation that results in cognitive impairment. Differential expression analysis is a popular technique for exploring the molecular underpinnings of neural aging, but technical drawbacks of the methodology often obscure larger expression patterns. Co-expression analysis offers a robust alternative that allows for identification of networks of genes and their putative central regulators. In an effort to expand upon previous work exploring neural aging in the marine model Aplysia californica, we used weighted gene correlation network analysis to identify co-expression networks in a targeted set of aging sensory neurons in these animals. We identified twelve modules, six of which were strongly positively or negatively associated with aging. Kyoto Encyclopedia of Genes analysis and investigation of central module transcripts identified signatures of metabolic impairment, increased reactive oxygen species, compromised proteostasis, disrupted signaling, and increased inflammation. Although modules with immune character were identified, there was no correlation between genes in Aplysia that increased in expression with aging and the orthologous genes in oyster displaying long-term increases in expression after a virus-like challenge. This suggests anti-viral response is not a driver of Aplysia sensory neuron aging.


Introduction
As an organ composed of long-lived cells, the brain is uniquely susceptible to the deleterious effects of aging, the outcome of which is often cognitive impairment [1,2]. Common hallmarks of brain aging include impaired metabolism, compromised proteostasis, mitochondrial dysfunction, and neuro-inflammation [3][4][5]. However, what causes these hallmark phenotypes is not well understood and still debated [6][7][8]. Due to the complexity of mammalian brains, invertebrate models are often employed for the study of aging in the nervous system at the neuronal level.
The marine model Aplysia californica is well suited for study of aging neurons. These mollusks live approximately one year in the wild and mariculture setting and have relatively simple nervous systems of only approximately 10,000 neurons grouped into well mapped circuits [9][10][11]. Studies on the life history and reflex behaviors of these animals in controlled environments have described weight loss and cognitive impairment in old age similar to that of other animals and humans [12][13][14][15][16]. Investigation of physiology and transcriptomics in easily identifiable neurons and neuron clusters during aging have elucidated neuronal correlates that underpin aging phenotypes at the behavior and organism level [17][18][19][20][21][22][23]. Indeed, it is the capacity for vertical integration of investigation, from molecular to behavioral, that makes Aplysia such an effective model for the aging nervous system. Previously, transcriptomic studies in aging sensory neurons of Aplysia identified signatures of common aging hallmarks, namely metabolic and proteostatic impairment [20,24]. However, due to the limitations of differential expression analysis (DEA), including log fold change thresholds and multiple comparison [25], the driving mechanism behind these transcriptional changes could not be identified. A common alternative analysis that can overcome these limitations is weighted gene correlation network analysis (WGCNA). WGCNA is able to lessen the impact of correction for multiple tests by comparing changes in groups of genes, called modules, as opposed to individual genes and does not employ the strict fold change thresholds used in DEA [25,26]. This method can capture network level changes that integrate the small, coordinated changes of many genes that would otherwise be below the thresholds employed in DEA [27]. Furthermore, WGCNA allows for the identification of putative central driver genes in co-expression networks and inference of putative function for genes that are undescribed or have yet to be experimentally verified via guilt-byassociation [26,28,29]. In this study we used WGCNA and eigengene network analysis to identify the central drivers of the transcriptional aging phenotype of Aplysia SN.

Experimental design
Sequencing data for this study were generated in our previous study [24]. RNA was extracted from A. californica Buccal S Cluster (BSC) and Pleural Ventral Caudal (PVC) sensory neurons across the adult aging spectrum (6-12 months).

Raw read processing
Raw, 150 base pair, paired end RNA reads from our previous study, available at the NCBI (PRJNA639857), were processed as described previously [24]. Briefly, raw reads were adapter trimmed and quality filtered using the BBDUK software [30] and then quantified by the Salmon software package [31] using the Aplysia californica reference transcriptome from the NCBI ftp site (AplCal3.0 GCF000002075.1).
Salmon derived transcript abundances were imported into the R statistical environment via the tximport R package [32]. Transcripts with a sum total abundance of less than 1 transcript per million (TPM) across all samples were considered not expressed and filtered out. The R package geneFilter was used to remove low variance transcripts using the function varFilter() with the default parameters var.func = IQR and var.cutoff = 0.5 to filter out transcripts with interquartile range (IQR) smaller than the median of all IQR in the expression data [33]. Surrogate variables were identified using the SVA R package [34]. Abundances were then variance stabilized using the vst() function from the DESeq2 R package, with surrogate variables incorporated into the model design and the blind parameter set to FALSE to account of said surrogate variables [35].

Co-expression analysis
Variance stabilized counts from BSC and PVC samples were used to construct sensory neuron type specific co-expression networks using the WGCNA R package. Briefly, Biweighted midcorrelation in WGCNA was used to construct adjacency matrices for each sensory neuron type independently. A soft power of 16 was used for both PVC and BSC networks to achieve at least 0.8 metric for scale free topology and to minimize mean connectivity (S1 Fig). Adjacency matrices were used to construct signed topological overlap matrices (TOM) for each sensory neuron type. The sensory neuron type-specific TOMs were centered and scaled and combined to make a consensus TOM by taking the minimum of the two modules. Transcripts were hierarchically clustered with the average method using a distance metric of 1-consensus TOM and a minimum module size of 30.
Initial module eigengenes were hierarchically clustered to determine module similarity. Modules with a branch height of 0.25 or less were deemed insufficiently different from their neighbors and merged. Module eigengenes for merged modules were then recalculated (S2 and S3 Figs). Consensus module eigengenes were then correlated to animal chronological age.
Individual transcripts in each module were then correlated to the module eigengene, referred to as module membership, as a measure of the transcript centrality.
Similarly, transcripts were correlated with chronological age, referred to as transcript-age significance (TAS), as a measure of the influence of chronological age on the expression of that transcript. Correlation of module membership with TAS of all transcripts within a module (MM-TAS) was calculated to describe the influence of chronological age on the module as a whole. Only modules with high magnitude (|Pearson cor| � 0.5) and high degree of significance (p � 0.01) of MM-TAS were considered for further downstream analysis. Module hub transcripts were identified using the WGCNA function chooseTopHubInEachModule().

Module enrichment analysis
Transcript sets of each module were tested for enrichment for Kyoto Encyclopedia of Genes and Genomes (KEGG) canonical pathways using the clusterProfiler R package [36].
All software used and version information is available in Supplementary Table S1 Table  and Supplementary Information S1 File. All scripts used for this analysis can be found in the following GitHub repository: [https://github.com/Nicholas-Kron/Kron_Cohort77_ CoExpression_Analysis]

Comparison to immune response in Crassostrea gigas
To assess whether the observed module immune signatures suggested mapping to KEGG orthology and the UNIPROT human proteome represented a bona fide molluscan immune signature, module transcript sets were compared to differential expression resulting from immune challenge in the pacific oyster Crassostrea gigas by Lafont et al (2020) [37].
The Aplysia RefSeq proteome for the current genome build (AplCal3.0, GCF_000002075.1) was downloaded from the RefSeq database. Similarly, the proteome for the C.gigas genome build used by Lafont et al (2020, oyster_v9, GCA_000297859.1) [37] was downloaded from the Ensemble FTP site. The Aplysia proteome was BLASTed against local BLAST database built from said C.gigas proteome using an e value cutoff of less than 0.001 and selecting only the top hit.
The resultant putative protein orthologs between Aplysia and C.gigas were then mapped to their respective gene and transcripts identifiers using the respective genome build gene feature format annotation files (AplCal3.0_genomic.gff version 1.21 and oyster_v9.49.gff3). The new Aplysia to C. gigas mapping file was used to determine the proportion of each co-expression module that mapped to genes differentially expressed in response to immune challenge either due to exposure to poly(I�C) priming, viral challenge, or both [37]. Enrichment for C. gigas orthologs in each module was calculated by Fisher's exact test with Bonferroni multiple test correction in R using the dhyper function. Modules were identified using the weighted gene correlation network analysis (WGCNA) R package. The most connected transcript, called the hub gene, is listed by its RefSeq identifier, as well as its BLASTx assigned human ortholog. Hub transcripts with a "-" in the Human Ortholog and Ortholog Name columns could not be mapped to any known human protein, and thus are of unknown function. metabolic failure such as Alzheimer's disease (ko05010), Huntington's disease (ko05016), and Parkinson's disease (ko05012). Orthologs also involved in anterograde or retrograde movement of cellular cargo featured prominently among the transcripts with highest module membership in this module. Metabolic enzymes associated with TCA cycle, glycolysis, and mitochondrial fatty acid beta oxidation were also prominent. Reactive oxygen species (ROS) detoxification enzymes and other mitochondrial homeostasis maintenance orthologs were also present ( Table 2).

Pink.
The pink module, which was positively correlated with age, was significantly enriched for KEGG pathways related to translation, such as ribosome biogenesis in eukaryotes (ko03008) and aminoactyl-tRNA biosynthesis (ko00970), and proteostatic mechanisms, such as lysosome (ko04142), protein processing in the endoplasmic reticulum (ko04141), and ubiquitin mediated proteolysis (ko04120, Fig 5).
Investigation of individual transcripts with the highest module membership for the pink module similarly revealed several transcripts annotated to human orthologs involved in transcription, translation, and ribosome biogenesis, as well as modulation of innate immunity and NFkB signaling (Table 3).

Orange.
The three KEGG pathways significantly enriched in the orange module, which was positively correlated with age, all participate in proteostasis, whether that is proper protein folding localized to the Endoplasmic Reticulum (ER) in the case of protein processing in the endoplasmic reticulum (ko04141) and N-glycan biosynthesis (ko00510), or in protein  degradation, in the case of lysosome (ko04142). Transcripts with highest module membership were associated with ER stress or the endoplasmic reticulum associated protein degradation (ERAD) pathway (Table 4).

Darkgreen.
For the darkgreen module, which exhibited an increasing eigengene expression trend until age 10 months after which the trend stabilized, KEGG enrichment analysis highlighted processes involved in nucleic acid metabolism, namely DNA replication (ko03030), Nucleotide excision repair (ko03420), and mismatch repair (ko03430) for DNA and RNA transport (ko03013) and RNA degradation (ko03018) for RNA (Fig 6).
Many of the transcripts with highest module membership were involved with DNA damage response, as well as formation of the nuclear pore, mRNA quality control and export, and immune signaling cascades such as the JAK/STAT cascade, NFkB signaling, and RIG-1 signaling (Table 5).

Greenyellow.
Pathways associated with viral infection and immune signaling dominated the KEGG enrichment results of the greenyellow module, which exhibited an increasing eigengene expression trend with age (Fig 7). Highly significant pathways included classical immune response associated pathways such as NF-kappa-Beta signaling pathway (ko04064), RIG-1-like receptor signaling pathway (ko04622), and NOD-like receptor signaling (ko04621). Transcripts with highest module membership in the greenyellow module mapped to human orthologs involved in the interferon and NFkB signaling pathways (Table 6).
A full list of KEGG enrichment results can be found in Supplementary Datasheet S1 Dataset. Full transcript sets and their respective module membership values can be found in Supplementary Datasheet S2 Dataset.

Module enrichment for C. gigas immune response genes
BLAST annotation of Aplysia proteins to the C. gigas proteome and subsequent annotation resulted in 22,715 unique Aplysia transcripts mapped to 10,112 unique C. gigas proteins. Of these 22,715 transcripts, 7,359 were present in one of the identified co-expression modules. Among the 1,547 genes marked as exhibiting coherent regulation profiles following priming with poly(I�C)and viral challenge in Lafont et al (2020) [37], 697 were also present in Aplysia co-expression modules. Enrichment tests for each module revealed that the greenyellow, darkgreen, and pink modules were significantly enriched for C. gigas immune response genes (Fisher's exact test, p < = 0.00385,  The RefSeq ID of each transcript is paired with a human protein ortholog gene symbol and name annotated by BLASTn mapping to the UNIPROT human proteome (UP000005640). The function of each ortholog is detailed in the final column. Correlation of transcript expression and chronological age, transcript-Age correlation (TAS), and significance (TAS p-value) are also listed. Transcripts are marked with " � " if they were among significantly downregulated transcripts in Kron et al 2020 [24] in one neuron type and with " �� " if in both. All subsequent tables are organized identically. Of the transcripts in the greenyellow, darkgreen, and pink modules, 17, 7, and 8 respectively exhibited a module membership greater than 0.8 for their respective modules and mapped to a C. gigas gene with a log 2 fold change of greater than or equal to 1 in response to immune priming and/or viral challenge in Lafont et al (2020) [37] ( Table 8). Of those, transcripts from the greenyellow and darkgreen modules mapped primarily to C. gigas genes with putative viral function. Furthermore, the greenyellow transcripts mapped to C. gigas genes assigned by   (2020) [37] primarily to the interferon-like and RIG-like receptor recognition pathways, while many darkgreen transcripts were assigned to the JAK/STAT signaling. Of those Crassostrea gigas transcripts that exhibited a greater than two-fold increase in expression 10 days after treatment with a viral analog, 84 had clear orthologs in Aplysia. Of those 84 orthologous transcripts in Aplysia, only two exhibited a significant increase in expression in aging sensory neurons in our previous study: an IRF8 ortholog and an uncharacterized protein (S4 Dataset).

Discussion
While enrichment analysis and eigengene expression profiles suggested strong overlap with our previous study, several modules and enrichment results identified many facets to the transcriptional dynamics in aging of these sensory neurons not detected in our previous DEA study [24]. Enrichment analysis and eigengene expression trend of the royalblue module strongly resembles that observed in Kron et al (2020) [24] of expression clusters with decreasing expression trends in aging. Key enzymes in glycolysis such as GPI, PGK1, PGAM2, GAPDH, and ENO1, as well as PDHB, which is part of the pyruvate dehydrogenase complex that links glycolysis to the TCA cycle, were among the transcripts with the highest module membership in the royalblue module. Key TCA enzymes DLST, IDH3G, and AOC2, as well as two members of the succinate dehydrogenase complex, SDHA and SDHC, which function as a nexus between the TCA and OXPHOS, were also present. Downregulation of these transcripts suggests decreased TCA cycle activity, and decreased NADH generation for use in OXPHOS. The presence of orthologs to OXPHOS complex assembly proteins SDHAF2, NDUFA2, and FMC1, as well as electron transport chain regulator ETRF1 among transcripts with high module membership may suggest further disruptions of OXPHOS at the complex level [41,71,73,75]. Mitochondria depleted of NADH exhibit impaired antioxidant capacity and increased ROS generation [258].
Several orthologs of major ROS detoxification enzymes, namely GPX4, SOD2, CAT, and PRDX5 are present in the downregulated royalblue module [40,54,55,80]. The further presence of PARK7, NXNL2, and ITCH, which stimulate antioxidant activity, in this downregulated See Table 2 for description of organization. Transcripts are marked with " � " if they were among significantly upregulated transcripts in Kron et al 2020 [24] in one neuron type and with " �� " if in both. Transcripts identified as orthologs to genes differentially expressed due to immune challenge in C. gigas are marked with a "!"in the first column.

PLOS ONE
Co-expression in aging Aplysia californica sensory neuron Nuclear pore formation and protein import into nucleus [194,195] (Continued ) a subunit of fast-inactivating A-type K channels, respectively, that repolarize neurons during firing and thus play important roles in membrane excitability, further suggests impaired neuronal function [63,74]. A facet of this module that was not captured in the expression clusters is the presence of many downregulated transcripts involved in cellular cargo transport. Orthologs of several proteins involved in retrograde transport, including DCTN6, CCDC151, and EFCAB1, and anterograde transport orthologs like BLOC1S1 and KIF3A, suggest disruptions in communication from soma to synapse and vice versa [38,44,51,56,66,76,259]. Downregulation of STAU2, which plays a key role in transport of RNA from the cell body to dendrites for local translation, suggests disruption of transcription/translation events necessary for long term memory in these sensory neurons [77,260]. Other processes in cellular cargo transport, such as ER to See Table 2 for description of organization. Common functions include nuclear pore formation, DNA damage repair, immune and inflammation signaling. The expression trend of this module increased linearly until age 10 months after which it stabilized, suggesting increasing activity of these in early age and stable, heightened activity in old age.   See Table 2 for description of organization. Common functions include ubiquitination, interferon signaling (IFN), and inflammation. The roughly monotonic increasing trend of this module's eigengene throughout aging suggests consistently increased activation of these processes in aging. https://doi.org/10.1371/journal.pone.0252647.t006 While the pink module was identified as having an increasing eigengene expression trajectory and a largely a transcriptional and proteostatic character by KEGG enrichment analysis similar to several clusters in Kron et al (2020) [24], the transcripts with the highest module membership, including the hub gene NFKBIA, suggest inflammation plays a central role in this module. Many of these upregulated transcripts mapped to human genes known to be induced by NFkB, including CSTL, BIRC3, FTH1, CYLD, and the hub gene NFKBIA [261][262][263][264][265][266]. Furthermore, several of these orthologs activate or permit NFkB signaling, including CSTL, BIRC3, GM2A, and NAT10 [83,97, [267][268][269][270]. The pink module also contains several dampeners of the NFkB signaling cascade and innate immunity, such as NFKBIA, CYLD, PDE12, SIGIRR, RIOK1, RIOK3 JKAMP, and TNIP1, likely to maintain homeostatic control and prevent over-inflammation [118,122,126,[129][130][131][132]271,272]. While an upregulation of NFKBIA would seem to suggest a decrease in NFkB signaling, the degradation of NKFBIA is a key step in the release of NFkB, allowing translocation to the nucleus [272]. This may suggest that these neurons are upregulating NFKBIA to keep up with growing rates of NKFBIA degradation as a result of increased NFkB signaling [122].
The pink module also contains upregulated translation modulators. This includes genes that promote rRNA maturation and recruitment such as EFL1, NAT10, RIOK1, RIOK3, and  HEATR1 [96,97,110,125,127]. Although known for its role in ribosome biogenesis and translation efficiency [97, 273,274], NAT10 is also responsible for N 4 -acetylcytoside (N4A) driven INF and NFkB inflammatory signaling via HMGB1 and NLRP3 inflammasomes, perhaps playing a dual role in this module [270]. Chronic, low grade inflammation as a result of N4A accumulation due to consistent activity of NAT10 may contribute to aging in these neurons.  [275][276][277]. However, this cascade has also been shown to be part of the antiviral response, specifically to prevent translation of viral RNAs [278,279]. Considering the preponderance of inflammatory genes in the pink module and the methodology that these animals were fed an ad libitum diet, the antiviral function of this cascade is the more likely here. Interestingly, this cascade also modulates synaptic plasticity and memory by inhibiting CREB in the hippocampus. Because CREB is essential for synaptic plasticity, increased CREB inhibition as a result of increased activity of the GCN1-EIF4A cascade during an inflammatory response may have knock-on effects that inhibit synaptic plasticity in these neurons [280].
Like the pink module, the orange module is similar to expression profile clusters found in Kron et al (2020) 281,282]. Sphingolipids, particularly ceramides, play central roles in pro-inflammatory signaling and ER stress, suggesting this module also participates in pro-inflammatory signaling [283]. SPLITC2 in particular has been shown to be upregulated by NFkB, suggesting NFkB signaling also plays a role in the orange module [281].
Several transcripts in the orange module regulate NFkB, such as BCL3 and BIRC3, similar to the pink module [83,149,150]. In addition, upregulation of the orange module hub gene CREB3L3 itself suggests neuro-inflammation as a central component of this module due to the role of this ortholog in the acute inflammatory response [135]. The similarities in the pink and orange modules suggest that each represents a different element of a proteostatic response to inflammation. Interestingly, the darkgreen module has a monotonic increasing eigengene expression trend early in the aging process and then stabilizes when the orange and pink modules enter their monotonic phase.
Many upregulated transcripts in the darkgreen module are orthologous to human genes associated with DNA damage response, including the hub transcript RPA2, RPA3, NSMCE1, ZRANB3, and TREX2 suggesting mounting DNA damage with age in these neurons [160,184,187,201,205]. The presence of several upregulated transcripts orthologous to genes critical to the formation of the nuclear pore such as NUP214, NUP98, and NUP88 as well as the stability and export of mRNA such as TREX2, GLE1, and TENT4A, and transcription and  [220,284], viral RNA sensor and activator of the interferon response pathway DDX58/ RIG-1 [285], and genes that support differentiation and maturation of immune cells such as EPSTI1, RBBP9, SSUH2 and GIMAP1 [163,182,183,218]. Aplysia immune hemocyte aggregates are known to play a pivotal role in neuron injury-associated inflammation and perhaps are similarly involved in the inflammatory response suggested by the orthologs present in the pink, orange, and darkgreen modules [286][287][288][289]. Curiously, TIPARP and PARP14 function to suppress and activate different cytokine signaling cascades depending on context and could be counted among other previously listed groups [166,167,198]. Upregulation of orthologs to the aforementioned viral RNA sensor RIG-1, sensor of viral entry PLA2G16, antimicrobial peptide MPEG1, and blocker of viral DNA replication BANF1 may suggest a role for this module in the detection of and initial response to viral infection in these neurons [191,193,204]. Interestingly, this may be further supported by the orthologs HENMT1 and MOV10L1 which mediate the formation of PIWI interacting RNAs (piRNA), noncoding RNAs that have antiviral activities in the innate immune systems of insects and possibly all eukaryotes [175,185,290,291]. The final module of interest, the greenyellow module, exhibits a monotonic increase in eigengene expression during the aging process, in parallel with the darkgreen module from age 6-9 months, and then parallel to the pink and orange modules thereafter. Transcripts with highest module membership in the greenyellow module all represent orthologs to elements of the Interferon (IFN) Mediated response to RNA viruses. Several innate immune mobilized antiviral zinc finger protein orthologs are represented in this module, including PARP12, ZC3HAVIL, TUT4, and three orthologs of ZNFX1, one of which is the hub transcript [249,252,292,293]. Other transcripts are orthologs of interferon-stimulated genes (ISG) that stimulate or facilitate the expression of other ISG, including ZNFX1, MRE11, DTX3L, CMTR1, and IRF1 [222,226,229,244,248]. Among these ISG orthologs are viral dsRNA sensors ZNFX1 and DDX58/RIG1 [222,285]. Other ISG orthologs upregulated in this module are IFN effector genes that have specific antiviral action: SAMD9 triggers anti-viral granule formation, TUT4 uridylates viral RNAs thus tagging them for degradation, SMG1 restricts viral replication, and TRIM2 prevents viral internalization [227,232,241,251]. The upregulation of TLR1 and STAT5b orthologs in the greenyellow module, which cooperate with STRA6 and JAK2 from the darkgreen, capture elements of the pattern recognition signaling cascade needed to mobilize an immune response [246,256,257]. A proper immune response also requires modulators to prevent over-inflammation or to dampen ISG expression once the viral challenge has been dealt with to allow clearance of waste. Several such immune regulator orthologs are also upregulated in the greenyellow module, including IL17RD which exists as part of a larger family of proteins that regulate interleukin (IL) and Toll-like receptor (TLR) signaling, ZNF598, and ASCC3 [234,239,243]. Whereas the pink module contains translation modulators, several of these upregulated transcripts are also orthologs of proteins known to have epigenetic and/or transcription attenuating effects, including PARP12, ZNFX1, and H3.1, possibly to mute global transcription in favor of anti-viral transcriptional programs [223,249,255].
Several of these transcripts also participate in NFkB signaling, promoting inflammation as part of the immune response. These include IRF1 and IRF8 which activate the Type-I Interferon response, CASP8, ASCC3 which is indispensable for NFkB signaling, and PARP12 which activates NFkB [243,249,254,294]. Ubiquitination dynamics are crucial for IFN and NFkB signaling [295], and many of the included upregulated transcripts are orthologs of genes that act as E3 ubiquitin-protein ligases such as DTX3L, ZNF598, TRIM2, and HERC4; as well as the Deubiquitinating enzymes VCPIP1. The multitude of orthologs involved in initiating the NFkB signaling cascade in response to viral infection in the greenyellow module may suggest this module may act upstream of the NFkB stimulated pink and orange modules.
While the discussed results appear to present a convincing picture of the function of these co-expression networks, it is important to approach these inferences with a degree of caution. Genes marked as orthologous by statistical means, such as BLAST or ghostKOALA, do not always perform the same functions in the target species as they do in the annotated species. Very few of these transcripts discussed have been independently investigated in Aplysia and their functions in Aplysia neurons and associated cells are not known. While several orthologs are generally understood to perform conserved roles across genera and their inferred function can be reasonably assumed, such as those involved in glycolysis or the TCA cycle, others involved in more idiosyncratic systems such as the immune response may have a large degree of divergence, especially when comparing the complex immune system of vertebrates such as human to that of evolutionarily distant mollusks.
Although the function of the orthologs discussed were characterized in human, mouse, C. elegans, and/or Drosophila, similar expression patterns to those captured by the pink, orange, darkgreen, and greenyellow modules have also been documented in other mollusks. Both Cathepsin L and MPEG1 orthologs have been shown to take part in the innate immune response of two species of abalone [296][297][298]. Immune challenge in oysters with viral RNA also resulted in upregulation of common NFkB and JAK/STAT signaling components observed in darkgreen and greenyellow modules, including orthologs to ZNFX1, Sacsin, and MPEG1 along with various IAP, TRIM, caspases, and IRF proteins [37]. To assess to what degree the transcripts inferred to play a role in the neural immune and inflammatory response via annotation to human orthologs represent a bonafide immune response, we further compared our module transcript sets to the transcriptional signature of an acute immune response to immune priming and RNA virus exposure in another mollusk, the Pacific oyster Crassostrea gigas [37]. The pink, greenyellow, and darkgreen modules were significantly enriched for orthologs to putative C. gigas immune response genes, supporting the notion that these modules do indeed represent a component of the Aplysia immune response. Indeed, several DE genes in the C. gigas immune response were orthologs of transcripts with high module membership in the pink, darkgreen, and/or greenyellow modules, such as the aforementioned orthologs of HERK4, BIRC3, IRF1, IRF8, MPEG1, DDX58/RIG-I, and most importantly all three orthologs of ZNFX1, including the greenyellow hub gene. Furthermore, several more transcripts, although not identified as orthologs of genes demonstrated to be DE in Lafont et al (2020) [37], do still map to the same human ortholog as genes identified in Lafont et al (2020) [37], such as sacsin. However, when comparing the list of significantly differentially expressed genes as a result of poly(I�C) in C. gigas to Aplysia orthologous genes significantly upregulated in sensory neuron aging as reported in our previous study, only two orthologs are shared. This suggest that, while the greenyellow and darkreen modules may capture an immune transcriptional program, viral infection and associated immune response are not drivers of transcriptional change in Aplysia sensory neurons. Moreover, this lack of correlation with oyster genes chronically upregulated after exposure to a virus analog indicates the need for caution in interpreting changes in module expression. Even when specific modules that show an increase in expression with aging include transcript members with anti-inflammatory or antiviral functions, these specific genes may not individually exhibit age-related expression changes. On the other hand, a majority of transcripts with high module membership in the royalblue, pink, and orange modules were identified as differentially expressed in our previous study, suggesting inflammation may result in proteostatic stress and mitochondrial dysfunction rather than infection.
Several transcripts from the royalblue module interact with the pink and orange modules, primarily as inhibitors. MTFS2 specifically inhibits the activity of EIF2AK3 from the orange module, which functions to maintain mitochondrial morphology, Ca 2+ homeostasis, and limit ROS production [299]. ITCH, via its interaction with TXNIP, promotes an antioxidant state, prevents pro-inflammatory signaling through ROS and NFkB, and prevents TXNIP from inhibiting glycolysis [300,301]. Similarly, PARK7 also inhibits NFkB signaling and the astrocyte inflammatory response [302,303]. Finally, GPS2 also exhibits strong anti-inflammatory activity [304,305]. Conversely, PDE12 from the pink module is known to suppresses mitochondrial translation and contribute to respiratory incompetence [306]. Perhaps the eigengene expression perturbation at 9 months of age common to these expression modules may represent a transition from a state of healthy metabolic activity and antioxidant capacity exemplified by the transcripts of the royalblue module towards an aged state typified by inflammation and protein dyshomeostasis suggested by the transcript sets of the pink and orange modules. The role of ROS represents a potentially interesting link between these modules.
Decreased activity of mitochondrial homeostasis and metabolism mechanism suggested by the royalblue module suggest mitochondrial dysfunction with age, a well-known source of ROS [307]. Furthermore, downregulation of several antioxidant proteins in the royalblue module would suggest decreased capacity for ROS defense. Overexpression of SOD2 and GPX have been demonstrated to inhibit NFkB activation [308,309], thus the downregulation of these and other antioxidants in the royalblue module may in fact potentiate NFkB activation as suggested by the pink and orange modules via increased ROS levels. While the role of ROS in NFkB signaling is complex and cell type specific, prolonged oxidative stress has been shown to increased NFkB signaling and contribute to a pro-inflammatory state [310,311]. Upregulation of FTH by NFkB, as seen in the pink module, has also been suggested to be a major component of the NFkB derived antioxidant response against high H 2 O 2 levels during chronic inflammation and oxidative stress [312]. Furthermore, ER stress resulting from increases in misfolded proteins due to oxidative stress, signatures of which were identified in our previous study and are also suggested by the orange module, has been shown to activate NFkB as well. Specifically, the activity of PERK, which is among the top orthologs in the orange module, plays a crucial role in activation of NFkB during ER stress [313]. In addition to downregulation of antioxidants, co-temporal upregulation of ROS producers like the two DUOX1 orthologs in the pink module suggest an increased in ROS, [88,134].
Indeed, chronic, low-grade inflammation has been suggested to be both cause and consequence of mitochondrial dysfunction and resulting metabolic impairment in neuronal aging [314][315][316]. These data suggest that chronic inflammation contributes to an increasingly oxidative environment in these neurons, exacerbating mitochondrial dysfunction that is known to occur in aging.  Table. Software. All Software and respective versions used for RNA sequencing read quality control and quality assurance, mapping and abundance estimation, and downstream analysis. (DOCX)

S2 Table. Transcript set overlap between co-expression modules (modules) and transcript expression profile clusters (clusters) from Kron et al (2020) [24].
Each cell represents the number of transcripts shared between respective module (row) and cluster (column). The "n" column represents the number of transcripts in a given module, and the "n" row represents the number of transcripts in a given cluster. Values in the "Sum" columns are row sums, e.g. the total number of transcripts in each respective module also present in the clusters. The "% of module" columns represent percentage values calculated by dividing the "Sum" columns by the total number of transcripts in a cluster set (1106 for B clusters, and 1198 for P clusters). Values in the "Sum" row are column sums, e.g. the total number of transcripts in each respective cluster that are also present in the modules. The "cluster %" row represent percentage values calculated by dividing the "Sum" row by the total number of transcripts among all modules (10012).