Transcriptional response of Caenorhabditis elegans when exposed to Shigella flexneri

In recent years, researchers have begun to use Caenorhabditis elegans as a potential animal model to study Shigella pathogenesis. This study aims to further develop this model using RNA-sequencing to understand which pathways/cellular characteristics are affected and potentially cause death in Shigella-exposed worms. We identified 1631 differentially expressed genes in Shigella-exposed worms (6 h exposure). A number of these genes encode proteins involved in fatty-acid β-oxidation (FAO), antioxidant defense and autophagy. The down-regulation of acyl-CoA dehydrogenases would impede FAO, reducing the overall energy to combat Shigella in the worm's intestinal tract. This is potentially coupled with the production of reactive oxygen species (ROS) that may not be fully quenched by antioxidant defense proteins, leading to damaged cellular organelles in the worm's intestinal cells. These cells may undergo autophagy to remove the mounting damage, but may eventually undergo cell


Introduction
Shigellosis (bacillary dysentery) is a highly invasive infection of the human colon that is responsible for approximately 190 million cases and 66,000 deaths each year [1]. The causative agent of shigellosis, Shigella, comprises four species that are further subdivided into serotypes: S. dysenteriae (15 serotypes), S. flexneri (19 serotypes), S. boydii (20 serotypes) and S. sonnei (one serotype) [2]. Due to their insusceptibility to stomach acids, only as few as 10-100 bacterial cells need to enter the colon in order to initiate pathogenesis [3]. Based on the distribution and prevalence of Shigella, S. flexneri has been identified as the leading cause of shigellosis in developing countries [3,4]. The severe damage caused to the colonic epithelial cells by the invading bacteria presents symptoms ranging from watery diarrhea to severe dysentery (characterized by fever, intestinal cramps escalating to bloody diarrhea with mucopurulent faeces) [5].
Due to its limited host range (humans and non-human primates), a feasible intestinal animal model has not yet been identified to study Shigella pathogenesis, thereby impeding the development of effective preventive and therapeutic measures [6]. The alternative animal models (murine pulmonary model and guinea pig keratoconjunctivitis model) that are currently available tend to utilize mucosal surfaces other than the intestine as the site of infection, questioning their clinical relevance (due to the symptoms presented in these models not reflecting the ones found in infected humans and the discrepancy in site of infection) [6][7][8].
Out of a number of model organisms currently being used to study host-pathogen interactions, the soil-dwelling roundworm, Caenorhabditis elegans, appears to be an attractive candidate to study Shigella pathogenesis for a number of reasons [6]. The intestinal cells of C. elegans share morphological similarities with human intestinal cells: finger-like microvilli on the apical side of intestinal cells that are anchored into a cytoskeletal terminal web made of actin and intermediate filaments. Also, the human innate immune system shares many characteristics with the immune system present in C. elegans, suggesting that C. elegans may utilize similar response mechanisms to counter a bacterial infection as its human counterpart [9]. Finally, its transparent anatomical structure, availability of self-fertilising hermaphrodites, rapid regeneration time and the requirement of a range of bacterial virulence factors to induce pathogenesis make C. elegans an ideal intestinal model to study the pathogenesis of Shigella and other major human bacterial pathogens (Pseudomonas aeruginosa, Salmonella enterica, Staphylococcus aureus, Enterococcus faecalis, Serratia marcescens, etc.) [10,11].
Recent experiments carried out on C. elegans exposed to Shigella show promise that the nematode could potentially be used as an in vivo animal model; the death of C. elegans exposed to the bacterium occurs via an infection-like process and is lethal in the presence of live bacteria carrying intact virulence plasmids [6,12,13]. Recent work carried out in our laboratory has shown that S. flexneri accumulates in the intestinal lumen of the nematode, produces outer membrane vesicles and eventually invades the intestinal cells [6]. Also, through the use of twodimensional differential in-gel electrophoresis (2D-DIGE), a number of proteins that were differentially expressed during the infection process were identified. Some of these differentially expressed proteins appear to be involved in the infection process, whereas others are induced as a protective mechanism [6].
In this study, we used high-throughput RNA sequencing (RNA-seq) to identify and quantify the differentially expressed genes in S. flexneriinfected L4 stage C. elegans. The differentially expressed genes identified in this study provide insight into the mechanism(s) involved in shortening the lifespan of C. elegans when infected with S. flexneri. We believe that this knowledge will lay a solid foundation for researchers dedicated to the development of this new in vivo model for studying S. flexneri pathogenesis and other bacterial pathogens.

Strains
The C. elegans wild-type N2 Bristol strain was obtained from the Caenorhabditis Genetics Center (CGC). N2 was cultured on nematode growth medium (NGM) agar plates seeded with E. coli OP50 and maintained at 22°C [14].
The bacterial strain, E. coli OP50, is a standardized food source for C. elegans in laboratory conditions [15]. Therefore, worms fed on E. coli OP50 provide the transcriptome profile of a healthy worm population that can be used to identify any changes in gene expression when the worms are exposed to other pathogenic bacteria/material [16]. S. flexneri 3b -SFL1520 (International Centre for Diarrhoeal Diseases Research Bangladesh-ICDDRB) was used in this study. Unless stated otherwise, E. coli and S. flexneri serotype 3b cultures were grown in Luria-Bertani (LB) broth and incubated overnight at 37°C and 30°C, respectively [6]. NGM plates used to expose C. elegans were seeded with either E. coli OP50 or S. flexneri 3b and, to ensure the expression of the genes present in the virulence plasmid, were initially incubated at 37°C overnight and kept at room temperature before transferring the worms.

Exposure of C. elegans to bacterial strains
L4 stage N2 worms grown on E. coli OP50 were washed thrice with M9 buffer and treated with 200 μg/ml Gentamycin for 3 h [14]. The worms were washed again three times with M9 buffer to remove any residual antibiotic and approximately 20,000 worms were transferred to three S. flexneri 3b and three E. coli OP50 seeded NGM plates. The plates were incubated at 22°C for 6 h. This time point was selected to ensure that RNA would be collected prior to the worms maturing to the adult stage. Otherwise, the adult worms would lay eggs, which would have hatched to L1 larvae. The presence of eggs and L1 stage N2 in the media would have greatly biased the results obtained from the RNA-seq analysis. After 6 h, the worms were washed with M9 buffer three times to remove the bulk of the residual bacteria present in the six samples. The worms were resuspended in 1.5 ml supernatant and transferred to siliconized microcentrifuge tubes and centrifuged for 3 min at 240 × g. The supernatant was removed and replaced with 1 ml sterile Milli-Q water. The samples were centrifuged as before and the excess supernatant removed. The tubes were placed in liquid nitrogen for 1 min and stored at −80°C.

Preparation of RNA-seq samples
The worm pellets were resuspended in 1 ml TRIZOL reagent. They were thawed and frozen in liquid nitrogen a total of six times. Afterwards, the suspension was vortexed (30 s) and put on ice (30 s) six times. The worm suspensions were transferred to new RNase-free Eppendorf tubes and allowed to sit at room temperature for 5 min. After the addition of 200 μl chloroform (mixed by inversion for 15 s and allowed to undergo phase separation for 3 min at room temperature), the samples were centrifuged at 14,000 × g at 4°C for 15 min. The RNAcontaining upper phase was transferred to new RNase-free Eppendorf tubes. Total RNA extraction was carried out using the RNeasy Mini Kit (Qiagen) according to the instructions provided by the manufacturer. The quality and the integrity of the total RNA extracted was assessed using a Nanodrop 1000 spectrophotometer (Thermo Scientific) and an RNA 6000 Nano Assay done using the Agilent 2100 Bioanalyzer (Agilent Technologies). Samples with an RNA integrity number (RIN) of 8 or higher were used in further analysis. 3 μg of total RNA at a concentration of 70 ng/μl (minimum volume 40 μl) was sent to the Ramaciotti Centre for Genomics, University of New South Wales, Sydney, Australia, for RNA-seq. Sequencing libraries were prepared using the TruSeq Stranded mRNA Library Prep Kit (Illumina). Libraries were sequenced using 2 × 75 bp sequencing on the Illumina NextSeq 500 platform.
Total RNA samples prepared similarly were used to carry out realtime quantitative two-step RT-PCR (qRT-PCR).

RNA-seq analysis
The quality of the 75 bp paired-end reads generated using the Illumina NextSeq 500 was assessed using FASTQC v. 0.11.3 [17]. The reads were first aligned to the E. coli K-12 and S. flexneri genomes to assess contamination, then aligned to the WBcel235 C. elegans genome assembly using HISAT2 under default parameters [18][19][20] (Supplementary file 1). Transcript abundances were quantified using feature-Counts v. 1.4.6 with strand-specific read counting [21]. The data normalization and differential gene expression analysis were performed using the limma+voom (using the voomWithQualityWeights function) in R v. 3.3.0 with a cutoff of false discovery rate adjusted P-value (FDR) < 5% and fold change (FC) > 2 (Supplementary file 2) [22][23][24]. Gene Set Enrichment analysis (GSEA) for Gene Ontology, KEGG and Reactome enrichments of the gene expression changes were carried out using the mroast function available in the limma package. The mroast function is utilized to identify whether processes with specified sets of genes are differentially expressed, while taking into account the directional annotation of the genes in the given process and the fold change of each said gene [22]. The enriched pathways identified via KEGG and Reactome pathway enrichment analyses required > 10 genes and FDR < 0.05 [25,26]. The heat map was generated using pheatmap v. 1.0.8.

qRT-PCR of C. elegans genes
qRT-PCR of C. elegans genes was carried out as previously described [11]. Briefly, total RNA from each sample was treated with TURBO™ DNase (ThermoFisher Scientific) according to the manufacturer's instructions. 500 ng of the DNA-free RNA was used in cDNA synthesis carried out using random hexamers (Roche) and Superscript™ III reverse transcriptase (ThermoFisher Scientific) at 50°C for 50 min. The synthesized cDNA was diluted 1:10. 1 μl of the diluted cDNA was used to carry out qPCR with POWER SYBR® Green PCR Master mix (Applied Biosystems) on a QuantStudio 3 Real-Time PCR System (Applied Biosystems) following the manufacturer's protocol except that the final reaction volume was set to 10 μl (genes tested from each biological replicate had three technical repeats). Each gene of interest was normalized to an internal control gene (cdc-42) and expressed as a log2 fold change of S. flexneri 3b infected samples compared with E. coli OP50 exposed samples [27].

Identification of differentially expressed genes in C. elegans exposed to S. flexneri
In order to identify the changes to the transcriptional profile of L4 stage N2 C. elegans when exposed to S. flexneri, we compared the gene expression change obtained from both infected and non-infected worms. The RNA samples utilized during this experiment were extracted from L4 stage N2 C. elegans exposed to either E. coli or S. flexneri for 6 h. With FDR < 5% and FC > 2-fold, we identified 1631 significantly differentially expressed genes (DEGs). Among them, 1119 and 512 are up and down-regulated genes, respectively ( Table 1, Fig. 1,  Supplementary file 3).
GSEA for Gene Ontology indicates that the expression of genes encoding a number of vital molecular functions was affected in the worms exposed to S. flexneri We next primarily focused our analyses on genes involved in processes that are known to have a significant impact on the worm physiology during a pathogenic bacterial infection [28,29].
3.2. Up-regulation of antioxidant genes in C. elegans exposed to S. flexneri Based on the GSEA for Gene Ontology, we observed that GO: 0016209 (antioxidant activity) was enriched ( Fig. 2A). Within the upregulated DEGs identified in Shigella-infected C. elegans, we found a number of genes encoding antioxidant defense proteins. This primarily Fig. 1. Volcano plot and heatmap of DEGs in C. elegans exposed to S. flexneri for 6 h. A: Volcano plot. Each colored point in the volcano plot represents a single gene. Plotted along the x-axis is the log2 fold change of each gene (comparison of the expression of each gene in worms exposed to S. flexneri with that of worms exposed to E. coli OP50). The y-axis represents the negative logarithm of the corresponding P-value of that gene. Up-(1119) and down-(512) regulated genes are shown in red and green (FDR < 0.05 and FC > 2), whereas genes showing no differential expression are shown in grey. B: Heatmap. The scale bar shows the z-score for a differentially expressed gene. N2_E1-3 designate RNA sample replicates obtained from C. elegans fed E. coli OP50 for 6 h, whereas N2_S1-3 are RNA sample replicates obtained from C. elegans fed S. flexneri for the same time duration. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 1 The number of DEGs in C. elegans when exposed to S. flexneri for 6 h. included a superoxide dismutase (sod) and three catalase (ctl) genes ( Table 2, Fig. 2B and Supplementary file 3). This list also includes a number of glutathione-s-transferases (gst), a glutathione peroxidase (gpx) and a thioredoxin (trx) gene, (Table 2 and Supplementary file 3) [30].
We were also able to identify that the gene encoding nicotinamide nucleotide transhydrogenase (nnt-1) was highly up-regulated in the infected C. elegans (FC: 13.48, FDR: 6.90 × 10 −7 ). This integral protein of the inner mitochondrial membrane uses energy from the mitochondrial proton gradient to generate large quantities of reduced nicotinamide adenine dinucleotide phosphate (NADPH). The NADPH generated is vital in the regeneration of reduced glutathione (GSH) from the oxidized glutathione (GSSG), to maintain a high GSH/GSSG ratio important for the reactive oxygen species (ROS) detoxification reactions carried out by GPX proteins; any C. elegans mutants lacking nnt-1 show higher susceptibility to oxidative stress generated via ROS [31,32]. Conversely, we were also able to identify a few gst genes to be downregulated (Table 2 and Supplementary file 3).
To confirm the differential expression of these antioxidant genes, we performed qRT-PCR on total RNA extracted from an independent S. flexneri infection to evaluate the transcript levels of three genes (sod-3, ctl-2, gpx-6 and gst-44). These four genes yielded similar results by qRT-PCR to those obtained by RNA-seq ( Fig. 3 and Supplementary file 5).
The up-regulation of such a high number of antioxidant genes within the S. flexneri-infected C. elegans suggests that there may be a rapid accumulation of various ROS within the host cells, after just 6 h of exposure to the pathogen. . These genes encode proteins known to be involved in the assembly of the autophagosome during autophagy [33,34].

Up-regulation of genes involved in autophagy
Among the up-regulated genes, we also identified a number of genes that encode proteins that participate actively in autophagy (Table 3). These genes predominantly encode aspartyl proteases (asp), cysteine proteases (cpr) and an asparaginyl peptidase, with some showing extremely high fold changes (Supplementary file 3).
To confirm the up-regulation of genes involved in autophagy, we performed qRT-PCR on lgg-1 (a candidate gene used as a marker to detect autophagy) [33,34]. lgg-1 gave similar results in both the qRT-PCR and RNA-seq analyses ( Fig. 3 and Supplementary file 5).
The potential differential expression of these proteins involved in autophagy, after only 6 h of contact with S. flexneri, suggests that certain cells, most likely the intestinal cells (the major contact between host cells and bacteria) in the infected C. elegans, may be starting to undergo degradation.

Differential expression of genes involved in fatty acid metabolism
Not surprisingly, certain pathways involved in the generation of energy (adenosine triphosphate [ATP]) to carry out essential cellular processes in C. elegans also appear to be affected by the presence of  Fig. 2A & D). However, expression of the genes encoding the acyl-CoA dehydrogenase proteins, acdh-1 (FC: 0.08, FDR: 6.32 × 10 −8 ) and acdh-2 (FC: 0.13, FDR: 3.45 × 10 −3 ), appear to be down-regulated in these worms.

Discussion
This study was carried out in order to probe the mechanism(s) that may play a key role in the death of C. elegans when exposed to S. flexneri. A large number (1631) of genes were identified as differentially expressed in the S. flexneri infected worms, even at an early time point of the infection (Table 1 and Supplementary file 3-4), suggesting an acute change in the overall transcriptional profile of the worm (Supplementary file 4). The genes with affected expression levels are involved in the immune response and the metabolic pathways that supply the energy required to combat the threat and whose overall  Fig. 3. Confirmation of the RNA-seq results using qRT-PCR. Expression of certain N2 differentially expressed genes was confirmed using qRT-PCR. All qRT-PCR results were normalized to the internal standard, cdc-42, before calculating the log2 fold change. All log2 fold changes in the expression of these genes obtained via RT-qPCR have a P < .05 (denoted by *). The error bars represent standard error of the mean. a The differentially expressed antioxidant genes were identified using limma (FDR < 0.05 and FC > 2). N/A: no DEGs were identified to be either up-or down-regulated for that antioxidant gene category. Table 3 DEGs involved in autophagy/cell death in C. elegans exposed to S. flexneri for 6 h. expression may possibly be 'hijacked' by effector proteins generated by the bacteria. It is feasible to hypothesize that this may cause the immune response generated against the bacteria to be less potent, and/ −or cause nutrients to be more readily available to the bacteria to carry out their own metabolic processes. We observed that expression of genes encoding acyl-CoA dehydrogenases were down-regulated. These proteins are essential for the initiation of fatty-acid β-oxidation (FAO) in the mitochondria (based on GSEA, the Gene Ontology term GO: 0006635 [fatty acid beta-oxidation] is significantly affected, as a higher proportion of genes in the geneset are down-regulated in S. flexneri exposed C. elegans) ( Fig. 2A & D) [35,36]. Acyl-CoA dehydrogenases catalyze the flavin adenine dinucleotide (FAD)-dependent oxidation of acyl-CoA to generate trans-2enoyl-CoA. The double bond in trans-2-enoyl-CoA is hydrated (via mitochondrial trifunctional protein [MTP] or crotonase) and the resulting L-3-hydroxyl-acyl-CoA is oxidized to 3-keto-acyl-CoA (via MTP or medium or short chain hydroxyacyl-CoA dehydrogenase). The 3keto-acyl-CoA undergoes thiolytic cleavage (via MTP or medium chain 3-ketoacyl-CoA thiolase) to produce a two-carbon chain-shortened acyl-CoA, one acetyl-CoA, one NADH and five reduced FAD (FADH 2 ), which are used during oxidative phosphorylation at the electron transport chain (ETC) to generate ATP [36,37]. The down-regulation of these genes, if reflected in reduced corresponding protein expression or activities, could result in a reduced rate of throughput of the FAO pathway. This could in turn hinder the capacity of the worm to generate sufficient energy to mount an effective immune response against the invading bacteria.
In contrast, the expression of a number of genes encoding lipases was highly up-regulated in the infected worms ( Fig. 2A & D and Supplementary file 3-4). Lipases are usually up-regulated during periods of starvation; they attack lipid droplets within the intestinal cells to break down triglycerides into fatty-acids, to be eventually used in FAO [33,38]. It may be that the infected worms up-regulate expression of lipase genes, but that the hypothesized slow rate of fatty acid oxidation, if confirmed, would allow the fatty-acids released by the lipase-activity to be utilized by the multiplying S. flexneri cells to synthesize phospholipids (an integral part of the bacterial membrane) and/or carry out their own FAO [39].
ROS are generated as a by-product of mitochondrial respiration. During high rates of aerobic respiration, electrons being transferred within the ETC may leak (mainly from ubiquinone, complexes I and III) into the matrix of the mitochondria, where they combine with oxygen (O 2 ) molecules to generate superoxide (O 2 -. ). Superoxide may dismutate to hydrogen peroxide (H 2 O 2 ), which in turn generates other oxidative species such as hydroxyl radicals ( . OH), peroxynitrite, maldonaldehyde and 4-hydroxynonenal (collectively referred to as ROS) [40]. The accumulation of high levels of ROS can cause severe damage to various cellular organelles by modifying DNA, proteins and lipids [40,41]. However, cells have evolved a number of antioxidant defense proteins to either detoxify the various ROS produced or repair the oxidative damage caused [40]. In C. elegans infected with S. flexneri, we observed the up-regulated expression of a number of genes encoding these antioxidant defense proteins ( Table 2, Fig. 2A-B and Supplementary file 3). If this is reflected in altered protein expression and function, the up-regulation of a large number of antioxidant genes suggests the presence of large quantities of ROS with potential to cause oxidative stress damage to the mitochondria and other cytoplasmic organelles.
Simultaneously, we observed differential expression of genes encoding proteins involved in autophagy. Autophagy, the bulk lysosomal degradation of cytoplasmic components, consists of a number of steps. These include the formation of a phagophore that encloses part of the cytoplasm (containing damaged organelles), through elongation and complete closure of its membrane to form an autophagosome, transport and fusion of the autophagosome with a lysosome, which initiates degradation of the enclosed cytoplasmic components [42]. Along with the up-regulated expression of a number of genes encoding proteins involved in the assembly of the autophagosome, we also identified a number of genes encoding cysteine and aspartyl proteases (cathepsins) ( Table 3, Fig. 2A & C and Supplementary file 3) [43]. These results could suggest that certain organelles within the cellular structures of the worm are extensively damaged (most likely due to the high accumulation of ROS) and are undergoing immediate degradation to preserve the viability of the cells. However, if the damage is extensive, or if the cells' ion (mainly Ca 2+ ) homeostasis is disrupted, these lysosomal proteases could also be released into the cytoplasm of the cell, leading to cell death [43]. Combined with the assault from multiplying bacteria, the lack of energy to stage an effective immune response, the accumulating ROS that damage cell (mainly intestinal) organelles, leading to cell death, these processes could be inferred as the cause of death in C. elegans exposed to S. flexneri.
A previous study used 2D-DIGE to identify the differentially expressed proteins in C. elegans exposed to S. flexneri for 24 h. qRT-PCR was performed on the genes that encode seven differentially expressed proteins that were identified using high stringency parameters: aco-1 (encodes an aconitase involved in maintaining iron homeostasis and is negatively regulated post-translationally by iron-levels); cct-2 (encodes a chaperonin that aids in the nuclear translocation of the transcription factor DAF-16. DAF-16 is the key transcription factor of the DAF-2/ DAF-16 insulin signaling pathway and is required for the expression of several antimicrobial genes); eef-2 (encodes the eukaryotic translation elongation factor that interacts with proteins involved in cell signaling, division, growth and development, ubiquitination mRNA splicing and energy metabolism); daf-19 (encodes a transcription factor involved in the p38 mitogen activated protein [MAPK] pathway and thereby involved in the innate immune response of the worm); hsp-60 (encodes a heat shock protein that regulates anti-bacterial immunity via the p38 MAPK pathway); unc-54 (encodes a protein involved in locomotion and pharyngeal pumping) and unc-41 (encodes a protein involved in locomotion) [6,[44][45][46][47][48][49]. It was shown that cct-2, daf-19, hsp-60 and unc-54 each exhibited a significant increase in transcript levels, whereas those of eef-2 and unc-41 were significantly decreased in transcript levels in the infected worms [6]. No statistical difference was observed for aco-1 [6].
In comparison, the RNA-seq results obtained in this current study, a 6-h infection, showed no significant changes in expression of aco-1, eef-2, daf-19, unc-54, unc-41 or cct-2. Only the expression of hsp-60 (FC: 0.39, FDR: 2.51 × 10 −6 ) was significantly down-regulated. It thus appears, not surprisingly, that the transcriptional, and translational, profile of C. elegans is likely to undergo further changes as the S. flexneri infection progresses, in order to regulate vital biological processes, e.g. maintaining iron homeostasis, mounting a stronger immune response. The worms would be responding robustly to neutralize the bacterial pathogens and their deployment of bacterial effectors to manipulate biological processes of the worm to facilitate leaching of nutrients that aid survival and multiplication of S. flexneri [6].

Conclusions
The results obtained from this RNA-seq analysis of C. elegans exposed to S. flexneri were utilized to investigate the possible cause(s) of death in the infected worms. We found that exposure to this bacterial species for a mere 6 h had a drastic impact on the overall transcriptional profile of the worm. Based on this transcriptional profile, we hypothesize that the presence of the bacteria could potentially lead to severe alterations to a number of metabolic pathways involved in generating ATP, thus reducing the ability of the worms to mount a strong defense. We also propose that there is a large accumulation of ROS within the cells of the worm that would cause extensive damage to the cell organelles, leading to autophagy and, finally, death of the cells. The overall damage to the cells, most likely the intestinal cells, may cause the eventual death of the worm. It is not yet clear which specific mechanism(s), if any, are utilized by Shigella to bring about such drastic changes within the worm. Therefore, we consider that our study opens new and exciting avenues of research to investigate such mechanisms to further improve our understanding of the C. elegans model of bacterial pathogenesis.