Wild-type Caenorhabditis elegans isolates exhibit distinct gene expression profiles in response to microbial infection

The soil-dwelling nematode Caenorhabditis elegans serves as a model system to study innate immunity against microbial pathogens. C. elegans have been collected from around the world, where they, presumably, adapted to regional microbial ecologies. Here we use survival assays and RNA-sequencing to better understand how two isolates from disparate climates respond to pathogenic bacteria. We found that, relative to N2 (originally isolated in Bristol, UK), CB4856 (isolated in Hawaii), was more susceptible to the Gram-positive microbe, Staphylococcus epidermidis, but equally susceptible to Staphylococcus aureus as well as two Gram-negative microbes, Providencia rettgeri and Pseudomonas aeruginosa. We performed transcriptome analysis of infected worms and found gene-expression profiles were considerably different in an isolate-specific and microbe-specific manner. We performed GO term analysis to categorize differential gene expression in response to S. epidermidis. In N2, genes that encoded detoxification enzymes and extracellular matrix proteins were significantly enriched, while in CB4856, genes that encoded detoxification enzymes, C-type lectins, and lipid metabolism proteins were enriched, suggesting they have different responses to S. epidermidis, despite being the same species. Overall, discerning gene expression signatures in an isolate by pathogen manner can help us to understand the different possibilities for the evolution of immune responses within organisms. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-08455-2.


Background
Antimicrobial resistance is a looming public health crisis as microbial pathogens evade drugs of last resort. When pathogenic microbes bypass the physical barriers of the host (e.g., an epidermis), an organism must rely on cellular and humoral defenses (e.g., the innate immune system) to protect itself from infections that can compromise organismal fitness. While all organisms will inevitably become infected over their lifespan, genetic variation in individuals will influence the outcomes of infections. Thus, our ability to link genetic variation to infection susceptibility may enable more sophisticated strategies for deploying antibiotics or help us enhance host immunity to treat microbial infection [1,2].
Genetic variation in organisms can create tradeoffs in susceptibility to different pathogens. The presence of non-synonymous single nucleotide polymorphisms (SNPs) in the coding sequence of immune effector genes has been shown to alter the propensity of individuals to become infected and develop disease [3]. For instance, a SNP in the TLR5 gene introduces a premature stop codon, creating a truncated form of the protein and impairing the ability of TLR5 to recognize the flagellin protein in flagellated bacteria [4]. Individuals carrying this SNP are more susceptible to infection by Legionella pneumophila, the vector for Legionnaires' disease, but not Salmonella enterica serotype Typhi, the microbial vector for typhoid fever [5,6].
Animal models have been critical to our understanding of the innate immune system. However, the reliance on canonical laboratory strains often ignores genetic variation in innate immunity and susceptibility to pathogens [7]. The soil-dwelling nematode Caenorhabditis elegans is a genetically tractable model organism to examine host-pathogen interactions [8,9]. C. elegans possesses a robust microbial defense system with evolutionarily conserved mechanisms in microbial pathogen detection, immune activation and pathogen killing (e.g., the induction of reactive oxygen species (ROS), deployment of antimicrobial and antifungal peptides, and expression of bacterial-binding C-type lectins) [10][11][12][13][14]. These evolutionarily conserved responses to infection make C. elegans an ideal model to study bacterial pathogenesis [9].
C. elegans wild isolates collected from geographically discrete regions are apt to vary in their ability to fight infections either because of a history of different selective pressures or simply genetic drift. Many discoveries regarding C. elegans host-microbe interactions have come from studying wild isolates, including identifying virus that infect C. elegans and the characterization of the C. elegans gut microbiome [15][16][17]. The genomic diversity of wild isolates has also provided a platform for quantitatively examining the genetic variation underlying complex traits, including behavioral differences in foraging as well as pathogen avoidance and susceptibility [18][19][20][21][22][23]. More recently,, studies have investigated variation in the transcriptomic responses of wild isolates to pathogens, such as Orsay virus [24].
In this study, we show that the C. elegans Hawaiian wild isolate CB4856 is significantly more susceptible to infection by the Gram-positive microbe Staphylococcus epidermidis than the N2 lab-conditioned strain. We used high-throughput RNA sequencing (RNA-seq) to identify differentially expressed genes in staged N2 and CB4856 worms exposed to pathogenic bacteria. The differentially expressed genes identified in pathogen-infected animals provide insight into isolate-specific transcriptomic responses of C. elegans to pathogenic bacteria. We believe that this knowledge will provide a foundation for better understanding the role of natural genetic variation in the C. elegans innate immune response to microbial pathogens.

CB4856 and N2 exhibit differences in susceptibility to S. epidermidis infection
During our studies on susceptibility to microbial pathogens we serendipitously cultured a bacterial contaminant that exhibited differences in a C. elegans isolate-specific manner. 16 S sequencing (see Materials and Methods) identified this bacterium as Staphylococcus epidermidis, an opportunistic pathogen and frequent cause of nosocomial infections [25]. We provided it with the strain designation EVL2000. Survival assays conducted with the EVL2000 S. epidermidis strain found CB4856 animals to be significantly shorter-lived compared to N2, with a median survival (LT 50 ) of 11.3 ± 1.6 days, whereas 50% of N2 animals were killed after an average of 18.8 ± 1.5 days (p < 0.001) (Fig. 1A, Table S1).
We then compared N2 and CB4856 susceptibility to other microbes that cause opportunistic infections, Staphylococcus aureus [26], Providencia rettgeri [27], and Pseudomonas aeruginosa [28]. In each case, the susceptibility of N2 and CB4856 to these pathogens was not significantly different (Fig. 1B-D, Table S1). We also raised animals on non-pathogenic E. coli OP50, the standardized food source for laboratory-bred C. elegans [29], to evaluate the possibility that the susceptibility of CB4856 to S. epidermidis was due to a general defect in lifespan. Interestingly, the lifespan of N2 worms was significantly shorter than CB4856, with 50% of animals dying at 12.0 ± 0.4 and 16.2 ± 0.6 days for N2 and CB4856 animals, respectively (p < 0.001) (Fig. 1E, Table S1). While there were differences in lifespan on E. coli, our data are consistent with the idea that the susceptibility of CB4856 to S. epidermidis is not due to a universal decrease in lifespan.
To determine if the differences in survival were unique to the EVL2000 strain, we infected N2 and CB4856 animals with an additional S. epidermidis strain (1191). CB4856 animals infected with the 1191 strain had a median survival of 5.7 ± 1.2 days, whereas 50% of N2 animals were killed after an average of 12.3 ± 0.9 days (p < 0.001) (Fig. 1F, Table S1). Lifespan of N2 animals infected with the 1191 strain was not significantly different from lifespan on E. coli OP50 (E. coli OP50 LT 50 = 12.0 ± 0.4 days, S. epidermidis 1191 LT 50 = 12.3 ± 0.9 days; p = 0.230). Together, our results indicate that CB4856, but not N2 is susceptible to S. epidermidis infection, and that this difference in susceptibility is observed across multiple strains of S. epidermidis. Since infection by the S. epidermidis EVL2000 strain provided the greatest difference in survival between N2 and CB4856, follow-up experiments used this bacterial strain.

Preference and avoidance behavior of N2 and CB4856 isolates depends on the microbe and time of observation
To determine if the difference in susceptibility could be explained in part by differences in animal behavior, we performed two-choice preference assays ( Fig. 2A). At the 1-hr time point there was a significant difference in S. epidermidis preference, with 100% of N2 and 69% of CB4586 animals on the OP50 lawn ( Fig. 2B) (p < 0.001). By 24 h the difference in preference was no longer significant and nearly 100% of worms from both isolates were found on the OP50 lawn. Preference assays in N2 and CB4856 animals given an option of E. coli OP50 and S. aureus, P. rettgeri or P. aeruginosa found differences in the initial preference behavior. At the 1-hr time point preference for S. aureus was not significantly different between the two isolates, with 81% and 68% of animals occupying the OP50 lawn for N2 and CB4856, respectively ( Fig. 2C) (p = 0.390). CB4856 animals were significantly more likely to prefer P. rettgeri at the 1-hr time point, with 52% of N2 and 19% of CB4856 animals on the OP50 lawn ( Fig. 2D) (p < 0.001). Preference toward P. aeruginosa trended in the opposite direction with 45% of N2 and 65% of CB4856 animals on the OP50 lawn ( Fig. 2E) (p = 0.051). Neither the N2 nor CB4856 genotype had significant changes in preference between the 1-hr and 24-hr time points for these three pathogens ( Fig. 2C-E).
We also examined avoidance behavior in N2 and CB4856 animals exposed to either E. coli OP05 or S. epidermidis (Fig. 3A). Animals of both genotypes exhibited little avoidance to E. coli. Both N2 and  (p < 0.001). Together, these results suggest that both genotypes prefer E. coli over S. epidermidis and both exhibit partial avoidance of S. epidermidis, albeit to varying degrees. Natural variation in chemosensation may partially explain these differences in behavior [19,30]. While the results are interesting and warrant follow-up, they likely have little impact on the difference in survival between N2 and CB4856 in response to monoculture on S. epidermidis.

Gene expression changes during infection cluster by genotype, rather than pathogen type
To identify gene expression changes due to microbial infection we performed RNA sequencing comparing CB4856 and N2 animals fed pathogenic bacteria versus non-pathogenic E. coli OP50 for 24 h. E. coli OP50 is a standardized food source for laboratory-bred C. elegans, and thus represents the nonpathogenic control diet. We analyzed gene expression in N2 and CB4856 animals exposed to either S. epidermidis, S. aureus, P. rettgeri or P. aeruginosa.
Although both C. elegans isolates showed similar susceptibility to three of these pathogens, we hypothesized geographical isolation might have resulted in variations in the immune response, which could be analyzed through transcriptional profiling. Principal component analysis (PCA) of gene expression found that N2 and CB4856 samples exposed to E. coli, the Gram-positive pathogen S. epidermidis or the pathogenic Gram-negative P. rettgeri clustered in an isolate-specific manner (Fig. 4A). N2 samples exposed to the virulent Gram-negative P. aeruginosa were also part of this N2-specific cluster, whereas CB4856 samples exposed to P. aeruginosa formed their own distinct cluster. Further, both genotypes exposed to the moderately virulent Gram-positive S. aureus formed distinct groups that did not cluster with other samples from the same genotype, nor did they cluster with each other (Fig. 4A).
We also performed PCA of gene expression using samples from a single genotype. For both C. elegans isolates S. epidermidis and P. rettgeri samples again clustered with E. coli, whereas. S. aureus samples formed a separate group along the first principal component axis (PC1). P. aeruginosa samples also formed individual clusters for   both isolates, with separation along the second principal component axis (PC2) in N2 and separation along both the PC1 and PC2 axes in CB4856 (Fig. 4B, C). Thus, it appears that both nematode genotype and microbial taxonomy contribute to the gene expression differences between N2 and CB4856.

Differential gene expression is largely unique in N2 and CB4856 animals after infection with pathogenic bacteria
We first sought to identify how gene expression changes in these two isolates, as a function of infection. Thus, we identified genes differentially expressed in response to all pathogens (for all analyses, differentially expressed genes (DEGs) were defined as those with an FDR-adjusted P-value < 0.05 and fold change ≥ 2, infected vs. control). We identified 273 and 117 genes that were differentially expressed in N2 and CB4856, respectively, as a function of infection. In either genotype the majority of DEGs were upregulated (N2, 236/273 = 86%; CB4856, 96/117 = 82%) (Tables 1 and 2, Table S2). Interestingly, only six out of the 384 total genes were differentially expressed in both N2 and CB4856 (N2-only, 267/273 = 98%; CB4856-only, 111/117 = 95%) (Fig. 5A, Table S2).
We examined the correlation of gene expression using the log 2 fold changes after infection. Expression levels were weakly, but positively, correlated between N2 and CB4856 (ρ = 0.097; p < 0.001) (Fig. 5B). A correlation of baseline gene expression was determined using normalized gene counts from worms exposed to the control bacteria (OP50 E. coli). In that context N2 and CB4856 had similar patterns of expression (ρ = 0.826), suggesting that gene expression in response to infection diverges from a common baseline. Taken together, these data suggest that while there is some conservation of expression in response to infection, the pattern is relatively weak between these two isolates.
We further examined the DEGs identified in both isolates. Of them, only Y65B4BR.1 has any functional annotation relating to infection. Previous data suggests that  RNAi knockdown increases susceptibility to infection by S. aureus [31]. A single gene, oac-14, was identified as differentially expressed in opposing directions, i.e., up in N2 and down in CB4856. RNAi knockdown delays or slows development, suggesting oac-14 contributes to organismal growth [32].
To take a more holistic approach, we used geneset enrichment analysis to analyze DEGs. The GO terms "Stress response: pathogen", "Transmembrane transport: lipid", and "Metabolism: lipid" were observed to be enriched in N2 (N2-only), whereas genes corresponding to the GO term "Metabolism: lipid" were overrepresented in CB4856 animals (CB4856-only) ( Table 3, Table  S3). Taken together, these results indicate that although there was little overlap in the transcriptomic response to microbial pathogen in N2 and CB4856 animals, both genotypes exhibit an enrichment of genes involved in lipid metabolism.
Finally, to examine the change in transcriptional networks, we cross referenced our DEG lists with a table of known C. elegans transcription factors (File S1). We found that in N2, the transcription factors cky-1, nhr-19, nhr-99, nhr-116, nhr-127, peb-1 and unc-42 were all upregulated in response to infection by pathogens, whereas in CB4856-infected animals, ham-1 expression was repressed, and klu-2 expression was induced (Table  S2). These results are consistent with the divergence of gene expression profiles between N2 and CB4856, suggesting that, upon infection, these isolates rely on different immune pathways.

N2 and CB4856 transcriptomic responses exhibit more overlap following infection with gram-positive than gram-negative pathogens
We then asked how gene expression profiles differed when we grouped pathogens by Gram-stain. We identified 3833 and 2565 genes that were differentially expressed in response to Gram-positive bacteria in N2 and CB4856, respectively. In either genotype the majority of DEGs were upregulated (N2, 2964/3833 = 77%; CB4856, 1962/2565 = 76%) (Tables 1 and 2, Table S4   CB4856 worms (Fig. 6A, Table S4). Further, the majority of DEGs identified in both genotypes were expressed in the same direction, while only 1% of genes were differentially expressed in opposing directions, i.e., down in N2 and up in CB4856 (14/1172 = 1%) (Fig. 6C, Table  S4). These results were reflected in a strong, positive correlation between N2 and CB4856 gene expression levels (ρ = 0.613; p < 0.001) (Fig. 6A). We again looked for transcription factors that were differentially expressed in both N2 and CB4856 as a function of infection by Gram-positive bacteria. In this list, there were 50 transcription factor genes that were significantly different from baseline (Table S4). As observed for the overall gene expression similarity, for each of the transcription factors, we observed a similar pattern for the direction of expression (e.g., genes upregulated in N2 were upregulated in CB4856, etc.).
These results highlight the dramatic differences when we compare the transcriptomic response of N2 and CB4856 worms infected with Gram-positive or   Gram-negative pathogens. Interestingly, exposure to Gram-positive pathogens resulted in differential gene expression profiles with considerable overlap between the two isolates, suggesting that there are conserved responses to some bacteria. In contrast, the transcriptomic responses of N2 and CB4856 animals to Gramnegative pathogens were relatively divergent, with few differentially expressed genes in common.
The majority of DEGs identified in both genotypes following either S. epidermidis or S. aureus infection were expressed in the same direction (S. epidermidis,   Table S6, Table S7). Lastly, we only identified nine transcription factors that were coordinately regulated in S. epidermidis-infected animals, whereas 346 transcription factors genes were differentially expressed in both genotypes in the same direction after S. aureus infection (Table S6, Table S7). Together, these results suggest N2 and CB4856 worms exhibit considerable similarity in their transcriptomic response to S. aureus, whereas there is much less overlap in their response to S. epidermidis infection, which may perhaps explain some of the difference in susceptibility to infection by S. epidermidis.
Taken together, these results indicate that the N2 and CB4856 transcriptomic responses to P. aeruginosa have more in common, whereas little overlap in differential gene expression exists when N2 and C4856 are infected with P. rettgeri. In fact, of all the gene expression profiles compared, P. rettgeri infection was the only one to have a negative correlation (ρ = -0.387; p < 0.001). These results are consistent with the responses to infection being driven both by the isolate and the microbe.

Geneset enrichment analysis identifies unique GO terms enriched in N2 and CB4856 worms exposed to S. epidermidis
Since N2 and CB4856 exhibited survival differences after S. epidermidis infection but not S. aureus infection, we were interested in identifying biological processes that may be enriched in these animals. We cross-referenced DEGs identified in N2 following S. epidermidis and S. aureus infection and identified 556 genes that were either exclusively expressed in response to S. epidermidis or were expressed in response to both pathogens but in opposing directions (e.g., upregulated in S. epidermidis and downregulated in S. aureus). The same analysis performed for DEGs in CB4856 animals identified 749 genes. Only 73 of these genes were differentially expressed in both N2 and CB4856 (N2-only, 483/556 = 87%; CB4856only, 676/749 = 90%) (Table S10).
We performed geneset enrichment analysis on genes that were exclusive to N2 (N2-only), genes that were exclusive to CB4856 (CB4856-only), and genes that were differentially expressed in both N2 and CB4856 animals (N2 & CB4856). Biological functions pertaining to the extracellular matrix, stress response, and metabolism were enriched in N2 and CB4856 worms following S. epidermidis infection ( Table 4, Table S11). The Gene Ontology (GO) terms "Extracellular material: collagen" and "Extracellular material: cuticlin" were enriched in N2-exclusive DEGs, whereas genes corresponding to the GO terms "Stress response: C-type Lectin" and "Metabolism: lipid" were overrepresented in DEGs exclusive to CB4856. Genes in the GO category "Stress response: detoxification" were overrepresented in both isolates, even if the individual genes were not necessarily the same. Genes that were differentially expressed in both isolates were enriched for a single GO term, "Stress response: pathogen". A separate group of genes that was expressed exclusively in CB4856 animals was also enriched for this GO term. Notably, we did not find significant enrichment of the GO term "Stress response: pathogen" in DEGs exclusive to N2 animals. Thus, when only considering genes that are differentially expressed in an S. epidermidis-specific manner, our results indicate a common enrichment of biological processes related to stress response and innate immunity in N2 and CB4856 animals. At the same time, we also observe an N2-specific enrichment of extracellular matrix genes and a CB4856-specific enrichment of lipid metabolism genes, suggesting a degree of divergence in the transcriptomic response of these isolates to S. epidermidis infection.

Comparing DEGs between N2 and CB4856 fed OP50 versus S. epidermidis
The genomes of N2 and CB4856 animals differ at approximately one polymorphism per kb and gene expression between these two isolates, including the expression of innate immune genes, is known to vary across similar developmental stages, even in the absence of infection [33][34][35]. Since principal component analysis found that genotype is the greatest predictor of isolate-specific differences in gene expression following S. epidermidis infection (Fig. 4), we performed an additional gene expression analysis to identify genes that were differentially expressed between N2 and CB4856 animals fed either the non-pathogenic E. coli OP50 or the pathogenic S. epidermidis. Using N2 expression as a baseline, we found 3384 and 1892 genes to be differentially expressed in CB4856 animals fed E. coli OP50 and S. epidermidis, respectively. Upregulated genes outnumbered those that were downregulated at a ratio of 2:1, regardless of whether the worms were fed OP50 or S. epidermidis (OP50 Upregulated, 2121/3384 = 63%; OP50 Downregulated, 1263/3384 = 37%) (S. epidermidis Upregulated, 1245/1892 = 66%; S. epidermidis Downregulated, 647/1892 = 34%) ( Table S12, Table S13).
Next, we cross-referenced genes that were differentially expressed between N2 and CB4856 animals in the OP50 and S. epidermidis conditions and classified them into four categories based on changes in gene expression, or lack thereof. Relative to N2, a total of 2280 genes were differentially expressed in CB4856 animals fed OP50 but were no longer differentially expressed when worms were fed S. epidermidis. This suggests that the final level of expression for these gene was similar in N2 and CB4856 animals following infection. A total of 1072 genes were differentially expressed between N2 and CB4856 animals fed OP50 and continued to be differentially expressed in the same direction in worms fed S. epidermidis, suggesting that they are likely uninvolved in the transcriptional response to the pathogen. A third group of genes, 788 in total, were exclusively differentially expressed between N2 and CB4856 animals fed S. epidermidis. Finally, we identified 32 genes that were differentially expressed between N2 and CB4856 animals fed either OP50 or S. epidermidis but in opposing directions, i.e., up in CB4856 fed OP50 and down in CB4856 fed S. epidermidis (Table  S14).
The isolate-and pathogen-specific expression pattern of genes from these latter two categories prompted us to perform geneset enrichment analysis to identify biological processes that may be differentially enriched in N2 and CB4856 animals following S. epidermidis infection. The 788 genes that were exclusively differentially expressed in CB4856 animals fed S. epidermidis were enriched for GO terms pertaining to immunity (Stress response: C-type lectin), cell signaling and transduction (Signaling: phosphatase; Signaling: tyrosine kinase), reproduction (Major sperm protein), posttranslational modification (Protein modification: glycoprotein), and cytoskeletal movement (Cytoskeleton: microtubule). From the 32 genes that were differentially expressed in opposing directions in CB4856 animals fed either OP50 or S. epidermidis geneset enrichment analysis identified GO terms corresponding to immunity (Stress response: pathogen) and the extracellular matrix (Extracellular material: collagen) ( Table 5, Table S15). Thus, by considering differences in gene expression between N2 and CB4856 animals fed non-pathogenic E. coli OP50, we have identified additional sets of genes that are enriched for GO terms pertaining to, among other categories, stress response, innate immunity, and the extracellular matrix.
Finally, to better understand the transcriptional regulators that might be significantly different upon infection by S. epidermidis, we looked for transcription factors that were differentially expressed only upon infection. There were 22 protein coding genes (1 pseudogene), that fit these criteria. Interestingly, 13 of the 22 were members of the nuclear hormone receptor (nhr) family. There are 266 nhr genes in our list of 1077 transcription factors. The incidence of nhr genes in our DEG list was significantly enriched (P = 0.007, Fisher's exact test, odds ratio 4.4).

Discussion
We conducted survival, behavioral and transcriptional assays using the C. elegans strains N2 and CB4856 to identify differences in pathogen susceptibility between two wild-type isolates. CB4856 was susceptible to infection by the Gram-positive pathogen S. epidermidis¸ whereas N2 was not. CB4856 and N2 were equally susceptible to a related species, S. aureus, and two non-Staphylococcus pathogens, P. aeruginosa and P. rettgeri. We confirmed the CB4856 susceptibility using a second isolate of S. epidermidis. What might underlie the differences in susceptibility? Initially we hypothesized that, perhaps, the N2 strain lacked the ability to detect S. epidermidis as pathogenic. However, both N2 and CB4856 prefer E. coli over S. epidermidis and both exhibit partial avoidance of S. epidermidis, albeit to varying degrees, suggesting they were averse to the bacteria. Further, we found robust changes in gene expression in response to the S. epidermis infection. Unexpectedly, we found that, despite P. rettgeri being equally pathogenic to both isolates, CB4856 animals preferred P. rettgeri to E. coli, while N2 did not. This observation was also highlighted by the fact that changes in gene expression were negatively correlated in these two isolates in response to P. rettgeri.

Isolates exhibit distinct transcriptomic responses to infection
The microbes used in this study are diverse both in terms of their genome size and known virulence mechanisms [36][37][38][39]. Principal component analysis of the sequencing dataset found that gene expression changes during infection clustered largely by genotype, rather than microbe, even in instances where susceptibility was equivalent. At a glance, this suggests that the strategies these strains have adopted for dealing with microbial infection have diverged since they were separated from one another.
When we looked for genes that were differentially expressed in response to infection, i.e., was identified in response to all microbes, there were hundreds, compared to the lists for individual microbes (thousands). Further, we found only six DEGs in common between N2 and CB4856, with those genes largely uncharacterized functionally. Overall, these results suggest that N2 and CB4856 respond to infection by altering gene transcription, but the core infection response represents only a small fraction of the total genes regulated, and that the core transcriptional response in each background has diverged. It also suggests there are still a number of genes involved in immunity that need to be characterized.
By categorizing the genes that were differentially expressed in response to all pathogens, we found an overrepresentation of genes involved in lipid metabolism in both N2-and CB4856-infected animals. Lipid metabolism has been shown to be required for proper innate immune system activation and pathogen defense [40]. These genes could play similar roles in both isolates and thus, while they may reflect a divergence in specific gene activation, they could represent a similarity in how N2 and CB4856 animals utilize lipid metabolism enzymes to facilitate the innate immune response.

Gram-stain level analysis suggests an absence of conserved master regulatory pathways
In many organisms, the Toll signaling pathway is activated by infection with Gram-positive bacteria and/or fungi, while the IMD pathway is activated in response to Gram-negative bacteria. Thus, we examined the results grouping microbes by Gram-stain, and then compared between isolates. Our results found that the majority of DEGs were strain-specific, consistent with a divergence of the responses to pathogens, but that, generally, there was a stronger conservation of responses to Gram-positive bacteria, compared to Gram-negative. Our results are consistent with the evidence that Toll and IMD signaling is diminished or absent in C. elegans as master regulators of immunity and suggests each isolate deals with Gram-negative and Gram-positive pathogens differently.
We examined the list of DEGs for transcription factors that could be contributing to the differential expression profiles in each of the isolates. There were hundreds of transcription factors regulated by infection with Grampositive pathogens, with 50 overlapping between the N2 and CB4856 isolates. In contrast, for Gram-negative pathogens, N2 had only nine differentially expressed transcription factors, while CB4856 had 97, with none in common.
One possible explanation for the large overlap in N2 and CB4856 transcriptomic responses to Gram-positive bacteria is the identity of the microbial pathogens. S. aureus and S. epidermidis are both members of the Staphylococcus genus, share a number of virulence factors and invoke overlapping host immune responses following infection [37,41,42]. Thus, the greater proportion of overlapping DEGs in N2 and CB4856 animals infected with Gram-positive bacteria may reflect a common response to Staphylococcus bacteria. Similar results for infection of N2 animals by Enterococcus faecalis and Enterococcus faecium suggest that, within species of pathogens, there are core transcriptional similarities [43].

Individual pathogens induce different responses in C. elegans isolates
Since N2 and CB4856 worms exhibited similar susceptibility to infection by P. aeruginosa, S. aureus or P. rettgeri, we posited that RNA sequencing would reveal considerable overlap in gene expression between the two isolates. Indeed, after infecting N2 and CB4856 worms with either P. aeruginosa or S. aureus we identified thousands of differentially expressed genes, a majority of which were coordinately regulated in both isolates, supporting the idea that infection by either of these pathogens invokes similar kinds of transcriptomic responses in N2 and CB4856 animals. However, animals infected with P. rettgeri exhibited very little overlap in gene expression with fewer than 20 genes differentially expressed in both isolates and expressed in the same direction.
It is possible that N2 and CB4856 worms have adopted two very different immune responses to P. rettgeri infection, with neither providing a perceivable survival advantage. But why do we not observe a similar lack of overlap in gene expression after infection with the Gram-negative pathogen, P. aeruginosa, which also showed similar lethality in N2 and CB4856 animals? Analysis of gene expression at regular intervals throughout an infection has demonstrated that the host immune response is a dynamic process which needs to be initiated and sustained to fight an infection [44,45]. At 24 h post-infection, the time point at which RNA was collected for gene expression analysis, animals infected with P. aeruginosa were much closer to dying than animals infected with P. rettgeri (N2 P. aeruginosa LT 50 = 3.8 ± 0.3 days, CB4856 P. aeruginosa LT 50 = 3.5 ± 0.3 days; N2 P. rettgeri LT 50 = 7.8 ± 0.7 days, CB4856 P. rettgeri LT 50 = 7.2 ± 0.6 days). Perhaps animals exposed to P. aeruginosa or P. rettgeri are at two different stages in combatting their respective infections and the overlap in gene expression between N2 and CB4856 animals grows larger as animals progress from earlier to later stages of infection.

Changes in gene expression by infection with S. epidermidis
We found differences in survival in response to S. epidermidis. There has been some limited investigation into how S. epidermidis kills C. elegans, however, in many cases, the research has used both S. epidermidis and S. aureus with relatively similar results, although there are reported strains of S. epidermidis that are not pathogenic [46]. Additional research suggests that S. epidermidis can accumulate in the intestine and relies in part on biofilm formation rather than stable secreted toxins to kill the host [47]. However, a global survey of transcriptional changes in response to S. epidermidis infection has not been performed. For each isolate we narrowed down the list of differentially expressed genes to those that were expressed in a S. epidermidis-specific manner and performed geneset enrichment analysis on sets of genes that were either exclusive to N2, exclusive to CB4856, or expressed in both N2 and CB4856 animals. We observed an enrichment of genes encoding extracellular matrix (ECM) proteins, such as collagens and cuticlins, exclusively in N2 animals. Cuticular collagens are an integral component of the C. elegans epidermal barrier and have established roles in the host epidermal immune response to pathogens [48][49][50]. Further, ECM proteins are likely involved in promoting longevity as changes in collagen gene expression have been used as in vivo biomarkers in C. elegans aging research [51]. While we cannot rule out the involvement of the host epidermidis in the different responses of N2 and CB4856 to S. epidermidis, we did not observe any epidermal abnormalities, such as blistering of the cuticle following infection. Nor did we observe swelling of the tail region, a phenomenon found in M. nematophilum infections [52]. Previous studies discovered S. epidermidis does not adhere to the cuticle and found that C. elegans mutants with altered cuticle structure exhibit wild-type sensitivity to S. epidermidis infection [47]. While our data suggest a role for collagen genes in isolate-specific responses to S. epidermidis, we do not believe they are solely responsible for the divergent responses to infection.
In DEGs exclusive to S. epidermidis-infected N2 animals as well as DEGs exclusive to S. epidermidis-infected CB4856 animals, our analysis detected an enrichment of genes encoding detoxification enzymes. These enzymes are primarily responsible for metabolizing xenobiotic substances and have demonstrated roles in the C. elegans immune response [53,54]. In both N2 and CB4856 animals, the DEGs were either members of the Cytochrome P450 (CYP) family or possessed UDP glucuronosyltransferase (UGT) activity. It is possible that these genes play similar roles in both isolates and thus may reflect a divergence in how N2 and CB4856 animals utilize detoxification enzymes to counteract the microbial infection.
Finally, we took advantage of the ability to look for genes with expression differences between N2 and CB4856 on control bacteria. This allowed us to sort differentially expressed genes in the S. epidermidis condition into those that started out different, but ended up being similar, and those that changed as a function of infection to end up at different levels of expression. A priori we might consider the former group as being less likely to underlie the differences in survival. In those analyses, approximately 90% of the genes were upregulated in response to infection.
Genomic variation is known to be an important contributor to infection outcomes in C. elegans. Allelic variation in npr-1 can influence survival, fecundity and avoidance of the microbial pathogens P. aeruginosa and N. parisii [21,22,63]. Similarly, genomic variation in cul-6 may be responsible for isolate-specific differences in susceptibility to Orsay virus [64]. It is as of yet unclear how much of the variation in survival that we observe is due to single genes or how allelic variation in genes like npr-1 or cul-6 impact the transcriptional differences we observe between N2 and CB4856.

Conclusions
In populations, genetic and genomic diversity drive phenotypes. Immune signaling is no exception, and yet, our understanding of C. elegans immunity has largely been driven by the analysis of the canonical laboratory strain, N2. We found that the N2 strain and CB4856 wild isolate were significantly different in their survival following S. epidermidis infection. It is important to note that S. epidermidis is unlikely to be a natural C. elegans pathogen and thus is unlikely to reflect real selective pressure in the wild. Nevertheless, it is possible that the relative abundance of related pathogens in their respective ecological niches resulted in the shift in immune capacity between the N2 and CB4856 lines. Unfortunately, without sufficient metagenomic sampling of the environment, this would be difficult to address. However, the robust differences observed suggests that isolating the animals and transferring them to the laboratory environment has fixed genetic variation that underlies microbe-specific susceptibility.
Perhaps more surprisingly, we found significantly divergent transcriptomic responses between these isolates, even when pathogen susceptibility was equivalent. These differences could come from protein coding variants in the molecules that detect pathogens or the transcription factors that coordinate immune defenses, non-coding variants in regulatory regions, or, more likely, a combination of the two. By narrowing down the regions of the genome that are important for the observed differences, we will be able to better understand how genetic variation contributes to different immune outcomes.

Identification of S. epidermidis strain
DNA was extracted from 1mL of fresh culture using a QiaQuick miniprep column (Qiagen) and the 16S ribosomal region was amplified using universal primers 27F (5'-AGA GTT TGA TCC TGG CTC AG-3') and 1492R (5'-CGG TTA CCT TGT TAC GAC TT-3'). PCR cycling conditions were as follows: (1) (3) 68 °C for 10 min. The amplified PCR product was purified, Sanger sequenced, and the NCBI BLAST tool was used to identify the species.

Survival assays
NGM plates were seeded with 40 µL of bacteria (E. coli, S. epidermidis, P. rettgeri, P. aeruginosa, or S. aureus) and incubated overnight at 37 °C. Plates were acclimated to room temperature and 30 late-stage L4 worms were added. Each day, living worms were transferred to fresh plates and the number of dead worms was recorded. Survival assays were performed in triplicate for each pathogen and worm strain at minimum.

Two-choice preference assays
NGM plates were seeded with a 10 µL dot of pathogenic bacteria (S. epidermidis EVL2000, P. rettgeri, S. aureus, or P. aeruginosa) on one half, a 10 µL dot of non-pathogenic E. coli OP50 on the other half and incubated overnight at 37 °C. 20 late-stage L4 worms were placed on blank plates for 30 min to remove any external bacteria while assay plates acclimated to room temperature. Worms were placed at the midline and their location on the plate was recorded at 1-hour, 2-hour, 4-hour, 8-hour, and 24-hour time points. To avoid ambiguity in the number of animals associated with the specific bacteria, only animals residing within the spot of bacteria were assigned "preference" for either the control (OP50) or pathogen.

Avoidance assays
NGM plates were seeded with a 20 µL dot of bacteria (E. coli OP50 or S. epidermidis EVL2000) and incubated overnight at 37 °C. 20 late-stage L4 worms were placed on blank plates for 30 min to remove any external bacteria while assay plates acclimated to room temperature. Worms were placed at the periphery of the bacteria and their location on the plate was recorded at 1-hour, 2-hour, 4-hour, 8-hour, and 24-hour time points. Animals that crawled onto the walls of the plate were excluded from analysis.

Exposure of nematodes to pathogen
NGM plates were seeded with 250 µL of bacteria incubated overnight at 37 °C. Approximately 2,000 L4 stage worms grown on E. coli OP50 were transferred to NGM plates seeded with bacteria and incubated at 20 °C for 24 h. After 24 h, the worms were washed with M9 buffer three times to remove any residual bacteria. The worms were resuspended in 100 µL of M9 buffer and mechanically disrupted in liquid nitrogen using a ceramic mortar and pestle. Frozen tissue was transferred to a microcentrifuge tube and 1 mL TRIZOL reagent was added. The tubes were flash frozen in liquid nitrogen and stored at -80 °C. Three biological replicates were collected for each worm and bacteria strain.

RNA isolation and sequencing
Frozen worm tissue was thawed, vortexed for 5 s and incubated at room temperature for 5 min. After the addition of 470 µL chloroform (mixed by inversion and phase separated for 2 min at room temperature), the samples were centrifuged at 15,000 RPM at 4 °C for 15 min. The upper aqueous phase containing RNA (approximately 600 µL) was transferred to a new RNase-free Eppendorf tube. Total RNA extraction was carried out using the Monarch RNA Cleanup Kit (New England Biolabs) according to the manufacturer's instructions. Total RNA was quantified using the Qubit (ThermoFisher Scientific) and RNA quality and integrity was assessed using the Agilent Tapestation 2200 (Agilent Technologies). Sequence libraries were prepared using the TruSeq Stranded mRNA Library Prep Kit (Illumina) and sequenced using single-end 1 × 75 bp sequencing on the Illumina NextSeq 550.

Statistical analysis
All statistical analysis was carried out using SigmaPlot 14.0 (Systat Software, Inc.). Log rank (Mantel-Cox) analysis was used for pairwise comparison of survival curves and to calculate median survival (LT 50 ). Chi-squared or Fisher's Exact Test was used to assess statistical significance of the bacterial preference and avoidance data. Spearman's correlation was used to determine the correlation coefficient between log 2 fold changes in gene expression. Statistically significant differences are defined in the figure legend and noted in the figure with asterisks.