Lipidome and Transcriptome Profiling of Pneumolysin Intoxication Identifies Networks Involved in Statin-Conferred Protection of Airway Epithelial Cells

Pneumonia remains one of the leading causes of death in both adults and children worldwide. Despite the adoption of a wide variety of therapeutics, the mortality from community-acquired pneumonia has remained relatively constant. Although viral and fungal acute airway infections can result in pneumonia, bacteria are the most common cause of community-acquired pneumonia, with Streptococcus pneumoniae isolated in nearly 50% of cases. Pneumolysin is a cholesterol-dependent cytolysin or pore-forming toxin produced by Streptococcus pneumonia and has been shown to play a critical role in bacterial pathogenesis. Airway epithelium is the initial site of many bacterial contacts and its barrier and mucosal immunity functions are central to infectious lung diseases. In our studies, we have shown that the prior exposure to statins confers significant resistance of airway epithelial cells to the cytotoxicity of pneumolysin. We decided to take this study one step further, assessing changes in both the transcriptome and lipidome of human airway epithelial cells exposed to toxin, statin or both. Our current work provides the first global view in human airway epithelial cells of both the transcriptome and the lipid interactions that result in cellular protection from pneumolysin.

secretes pore-forming proteins during infection that aid in bacterial colonization. This pore-forming toxin, named pneumolysin (PLY), is a member of the cholesterol-dependent cytolysins (CDC) family secreted predominately in gram-positive bacterial pathogens 7 . Pneumolysin is secreted as a soluble monomer, which upon binding to host cell membranes via cholesterol, can oligomerize to form β -barrel pores and enter cells 8,9 . CDC infection is quite complicated, resulting in responses not only at the cellular membrane but also at many organelles inside the cell. Once the membrane is damaged, the cell works to repair the pore, removing the lesion in order to maintain its viability 10 . Yet, eventually the damage overwhelms the cell, ultimately resulting in cell death.
Statins are competitive inhibitors of 3-hydroxy 3-methylglutaryl coenzyme A (HMG-CoA) reductase, a key enzyme regulating cholesterol biosynthesis 11 . Due to their ability to inhibit cholesterol production, we began to investigate their use in protection from cholesterol-dependent cytolysins. More than 20 independent clinical epidemiological studies have been completed, suggesting that statins have a strong beneficial effect against pneumonia and sepsis related mortality 12 . Yet, statins' beneficial role is still considered controversial and additional studies have suggested these so-called protective effects are due to misinterpretation of data, with one study actually showing an increase in infections in statin users 13 . In contrast, laboratory studies using animal models of infection have shown protection against bacterial infections under statin use, mainly focusing on Staphylococcus aureus [14][15][16] . We decided to test whether these potentially protective statins dependent effects could occur in the airway epithelium, the main physiological target of bacterial infections. Indeed, we have found that statins pretreatment led to protection in vitro from pneumolysin cytotoxicity in human airway epithelial cells. To further understand the global response networks behind this protection, we characterized both the transcriptome and lipidome of human airway epithelial cells exposed to toxin, statins or both. Due to the fact that we were exposing human airway cells to statins, competitive inhibitors of 3-hydroxy 3-methylglutaryl coenzyme A (HMG-CoA) reductase, a key enzyme regulating cholesterol biosynthesis that results in changes to cholesterol dynamics, it was logical to assess changes in lipid composition 11 . Currently, there are only a limited amount of reports involving the lipidomic analysis of statins, with a majority carried out with human plasma samples [17][18][19] . This study is the first to analyze lipid profiles of both statin-treated and pneumolysin-treated human airway epithelial cells.

Results
Statin pretreatment led to global changes in plasma membrane lipid composition. To help determine the molecular landscape characteristics of the affect of statins on pore-forming toxins, we decided to perform two types of global analysis on these lung epithelial cells. Since statins affect an acetyl-CoA driven pathway that plays a central role in lipid biosynthesis, we performed shotgun lipidomics on HBE1 cells grown under typical conditions, focusing on a rapid comparison of lipid profiles using nanoelectrospray ion trap tandem mass spectrometry (nanoESI-MS). We chose to use an immortalized cell line over primary lung epithelial cells due to the high technical noise present in this technique that would be hard to discern from biological variation. This method produced over 700 distinguishable ions per mass spectrum in positive mode and over 600 distinguishable ions per mass spectrum in negative mode. Although not all of these spectra could be identified due to the lack of liquid chromatography separation and subsequent lipid adduct ambiguity, 132 unique lipid structures plus 22 metabolites were found using an in-house LipidBlast database (Table ST1). Of these, 89 lipid structures were quantitated to determine altered lipids between the four treatment groups: EC (no treatment), EP (pore-forming toxin only), SC (statin only), and SP (both statin pretreatment and pore forming-toxin). This is shown in Table ST2 and Fig. 1. We found that 24 lipids were significantly changed in airway epithelial cells, using ANOVA and a P value cutoff of less than 0.01 (Table 1). Of the 89 quantifiable lipids, a majority of the changes were seen in major constituents of biological membranes: phosphatidic acids (PA), phosphatidylcholines (PC) and phosphatidylethanolamines (PE). Looking in more detail, PLY treatment led to a shift in lipid composition away from normal levels, concentrating on increasing saturated and shorter chained lipids (i.e. Addition of lysoPC 18:1 did not restore toxin dependent effects. C18:1 lysophosphatidylcholine or lysoPC(18:1) was increased during toxin treatment and statin pretreatment led to a statistically significant decrease, on par with normal levels of this lipid (Fig. 3). Lysophosphatidylcholines are derived from phosphatidylcholines by partial hydrolysis via phospholipase A2 enzymes 20 . Due to its commercial availability, we decided to investigate whether lysoPC(18:1) supplement can abrogate the protective effect of simvastatin. HBE1 cells were pretreated with simvastatin (100 nM) accompanied with various doses of lysoPC(18:1) for 24 hours, then treated with pneumolysin (400 ng/mL) for 4 hours. By ATP release assay, we found that lysoPC(18:1) supplement did not abrogate the simvastatin-mediated protection against pneumolysin (Fig. 3). Interestingly, lysoPC(18:1) treatment led to an increase in cell viability, although not statistically significant, in both the simvastatin supplemented and control cells. Additionally, we showed a similar affect in simvastatin treated cells using the MTS assay and this was included as supplemental figure S2.
Transcriptome analysis indicated a strong role for early growth response genes. In order to better understand the cellular mechanism behind statin cellular protection, we completed transcriptome analysis. Lung epithelial cells from three patients with no known lung diseases were used in this study and grown in proliferating conditions. Primary cells were used in this experiment due to the low inherent technical variation present in mRNA-seq. As before, 24 hour pretreatment with simvastatin or mock control was followed by four hours of pneumolysin exposure or mock control. This resulted in four different classes from three different patients: vehicle control, only statin, only toxin, and both statin and toxin. In order to gain a better insight into the inherent nature of these twelve libraries, unsupervised clustering was done by way of multidimensional scaling or MDS plot. One of the main issues when working with primary human tissue samples is the source of noise due to biological variation. This was exemplified in the MDS plot, shown in figure S1; the twelve libraries separated out into three clusters based on their human sample identification and not into four clusters based on their treatment class. This indicated an inherent issue with the data; the variation between patients in each treatment class is high, reducing the power of this experiment and making statistical testing more difficult. Therefore, in order to compensate for this nuisance factor, and due to the fact that this is a paired study design, we can compare the treatment to the control for each set of lung epithelial cells, effectively subtracting out the baseline differences between the patients. Using this additive model and fitting a generalized linear model within the edgeR package, we were able to determine statistically significant genes. At a P value of 1% or less, the none versus toxin comparison had 430 significant transcripts, the none versus statin comparison had 318 significant transcripts and the none versus both statin and toxin comparison had 2924 significant transcripts. The top twenty genes ranked based on fold change over untreated cells for each of these comparisons are shown in table 2 and qPCR verification is shown in figure S3.  Table 1. Statistically significant lipids from mass spectrometry data. Intensities were normalized by class. Statistically significant differences were assessed by one-way ANOVA with post hoc Dunnett's correction and represented by asterisks (ns P > 0.05, * P < 0.05, ** P < 0.01, *** P < 0.001).
Scientific RepoRts | 5:10624 | DOi: 10.1038/srep10624 As was expected, members of the HMGR pathway were pulled out in the statin alone treated group, such as farnesyl diphosphate synthase (FDPS), 3-hydroxy-3-methylglutaryl-CoA (HMG-CoA) synthase 1 (HMGCS1) and acetyl-CoA carboxylase beta (ACACB). Interestingly, there were many genes in common between the toxin and statin alone treatments. Surprising to us, was the increase of both TNF and IL-8 in the statin alone treated group, at levels similar to those seen in the toxin alone treated group, implying a similar inflammatory response to statin exposure. In the toxin alone treated group, many of the top hits were transcription factors, such as EGR1 (Early growth response gene 1), EGR2 (Early growth response gene 2) and ATF3 (Activating transcription factor 3). EGR1 and EGR2 are members of the immediate early gene family and are induced by a variety of different stimulants: growth factors, neurotransmitters, hormones, stress and injury. Stress response genes were also found in the toxin alone treated group. ATF3 is a member of the cAMP responsive element-binding (CREB) protein family of transcription factors and involved in the cellular response to stress. CREBRF (CREB3 regulatory factor), a gene also up regulated in the toxin alone treated cells, acts as a negative regulator of the unfolded Human bronchial cells (HBE1) underwent the following treatments, for a total of seven replicates: no treatment/vehicle control (EC), pneumolysin (400 ng/mL) alone for 4 hours (EP), simvastatin (1 μ M) 24 hour treatment alone (SC) and both simvastatin (1 μ M) 24 hour pretreatment and pneumolysin (400 ng/mL) for 4 hours (SP). Intensities were normalized by class. Mean ion intensities with standard errors (boxes) and quartile ranges (whiskers) are shown. Positive mode and negative mode refer to the setting on the mass spectrometry under which this data was obtained. Statistically significant differences were assessed by one-way ANOVA with post hoc Dunnett's correction and represented by asterisks as compared to Control.
protein response (UPR), an endoplasmic reticulum stress response 21 . Notably, interleukin 20, as well as interleukin 8, were up regulated in both the toxin alone treated group and the combined treated group. Interleukin 20 has been shown to promote S. aureus infection, another pore-forming toxin producing bacteria, by down regulating IL-1β -and IL-17A-dependent pathways 22 . Most of the genes found in the both statin and toxin treated group are a combination of the two individual treatments. For example, the fold change associated with both IL-8 and TNF are equal to the addition of the levels of the two Mean ion intensities with standard errors (boxes) and quartile ranges (whiskers) are shown. Positive mode and negative mode refer to the setting on the mass spectrometry under which this data was obtained. Statistically significant differences were assessed by one-way ANOVA with post hoc Dunnett's correction and represented by asterisks as compared to Control. B.) HBE1 cells were pretreated with simvastatin (100 nM) with or without indicated doses of lysoPC 18:1 for 24 hours and then challenged with pneumolysin (400 ng/mL) for 4 hours. Cellular ATP release assays were used to assess cell viability. Statistically significant differences between groups were assessed by one-way ANOVA with post hoc Tukey's correction and represented by asterisks. Error bars represent S.E.M. of 3 experiments. individual treatments. Of note is the presence of phospholipase A2, group III (PLA2G3) in the both treated group. Secreted phospholipase A2 members have a wide variety of biological functions, including generating free fatty acids and lysophospholipids, and play a strong protective role in host defense against bacteria and viruses 23 . Interestingly, many of the genes regulated by EGR1 and EGR2 are also downstream of SREBP1 (Sterol regulatory element-binding transcription factor 1) and SREBP2 (Sterol regulatory element-binding transcription factor 2). Sterol regulatory element-binding proteins are a family of transcription factors involved in both lipogenesis and cholesterol homeostasis by controlling the expression of various enzymes required for cholesterol, fatty acid, triacylglycerol and phospholipid synthesis 24 . SREBPs are directly regulated by cholesterol levels and upon sensing a decrease, SREBPs are cleaved from the endoplasmic reticulum and allowed to translocate into the nucleus 24 . Through searching for transcription factor binding sites within 500 bp of the transcriptional start site (TSS) of a list of up regulated genes (at a P value of 1% or less) from the Statin only treated group, 67 binding domain motifs, with 47 unique transcription factors, were found. Of these, SREBP1 and SREBP2 were pulled out at the top of the list. This was completed using a plugin for cytoscape called iRegulon 25 and the output is included in the supplemental as table ST4. We combined genes regulated by SREBP1/SREBP2 and EGR1/EGR2 as well as fold change expression values for all three treatment groups compared to untreated controls and visualized this data in fig. 4.

Gene Symbol
Gene Definition  Cluster analysis results reinforced the involvement of fatty acid metabolism. A common issue with genome wide RNA expression analysis involves interpreting the data in order to gain biological insights. Gene Set Enrichment Analysis (GSEA) helps in finding ways to solve this computational pigeonhole by providing an easy to use program that aids in finding overall themes of any large dataset 26 . Created by the Broad Institute, GSEA takes advantage of the knowledge of the scientific community and looks at the data in regards to sets of related genes (gene sets or GS). Due to this grouping into clusters, an increase in statistical power is gained and potentially changes the outcome of an otherwise unfruitful experiment.
Looking at the comparisons, many gene sets are pulled out as being significantly enriched. In both the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets, none versus statin, toxin or both statin and toxin comparisons, many of the same gene sets are pulled out in all three cases. All three cases have quite a large number of hits when using the 25% FDR cutoff recommended for hypothesis generation by the GSEA authors. In the GO work, statin alone group has 48 gene sets, toxin alone group has 53 gene sets and the both treatment group has 46 gene sets. In the KEGG work, statin alone group has 49 gene sets, toxin alone group has 43 gene sets and the both treatment group has 44 gene sets. The top ten hits for each comparison from this analysis are shown in table 3 and table 4. As expected, due to the changes seen in the lipidomics work, many pathways involving fatty acid metabolism were found in both GO and KEGG gene sets (i.e. GO gene set lipid metabolic process and KEGG gene set biosynthesis of unsaturated fatty acids). Interestingly, translation was exclusively found in the toxin GO analysis, implying a dependency on protein production. Similar to the individual gene analysis, the gene ontology analysis shows an inflammatory response in the statin treated group, as well as the toxin alone and the both treated groups. Ultimately, many of the same pathways are found in all three comparisons, making it hard to discern out the specific biological processes involved in each Figure 4. Pathway analysis showed a strong role for early response genes and sterol regulatory elements. Lung epithelial cells from three patients with no known lung diseases underwent the following treatments: no treatment/vehicle control, pneumolysin (400 ng/mL) alone for 4 hours, simvastatin (1 μ M) 24 hour treatment alone and both simvastatin (1 μ M) 24 hour pretreatment and pneumolysin (400 ng/mL) for 4 hours. Iregulon cytoscape plugin was used to interrogate the proximal promoter (500 bp within the TSS) of statistically significant transcripts (p value adjusted for false discovery by Benjamini-Hochberg cutoff of 0.1 or less) upregulated in response to treatment as compared to no treatment. EGR1 (Early growth response gene 1), EGR2 (Early growth response gene 2), SREBP1 (Sterol regulatory element-binding transcription factor 1) and SREBP2 (Sterol regulatory element-binding transcription factor 2) were all found to have putative binding sites to all transcripts shown (pink nodes). Normalized mRNA-seq expression levels, averaged across all three patients, are shown via a pink to red gradient for all transcripts as fold change values, comparing treatment indicated to no treatment/vehicle controls.
Scientific RepoRts | 5:10624 | DOi: 10.1038/srep10624 treatment. In order to assist with this task, we completed leading edge analysis within GSEA on all three comparisons in order to better understand the main genes at work in each case. This type of analysis focuses on genes that are present in a subset of gene sets (chosen by the user; we focused on gene sets with a FDR of 25% or less) and assists in pulling out key genes important in each comparison. Although again many of the same genes were pulled out in all three comparisons, there are a few notable differences. Various members of the phospholipase A2 (PLA2) superfamily were pulled out in the none versus statin comparison and in the none versus both comparison (PLA2G3, PLA2G6 and PLA2G10), as seen in the single gene analysis. Again, these are of particular importance since there is a considerable body of evidence supporting the antibacterial function of PLA2 family members. Although it was not a top gene hit, in the none versus both comparison, IL-6 was found in 6 gene sets. IL-6 plays an important role in  host defense against bacterial infection; using IL-6 -/-mice, it has been shown that IL-6 is required for resistance against S. pneumoniae 27 . Pore-forming toxins result in a wide variety of cellular responses and GSEA was able to pick out some important pathways involved.

Discussion
In the present study, we have undertaken the first of its kind to investigate the changes in both lipids and RNA expression in lung epithelial cells exposed to a pore-forming toxin. Based on our data showing a protective effect from statin-pretreatment (as seen in Fig. 3B), we included this treatment in our current study in order to better understand the complex cellular processes associated with statin treatment.  Statins were originally discovered by Dr. Akira Endo in the 1970s as a molecule produced from the fungus Penicillium citrinum as a means of self-preservation from other microbes and shown to inhibit the enzyme HMG-CoA reductase, leading to a decrease in cholesterol production 28 . It has been known for years that bacterial-fungal communities create environmental conditions that assist in controlling the growth of other microbes. Furthermore, due to the fact that the pore-forming toxin of interest in this study, pneumolysin, is a cholesterol dependent cytolysin, it is logical to believe that statin pretreatment may have an effect on PLY toxicity. There are many intrinsic cellular pathways that have been shown to provide protection to pore-forming toxin (PFT) cytotoxicity, including mitogen-activated protein kinase (MAPK), Unfolded Protein Response (UPR), autophagy and calcium-dependent endocytic pathways, although the molecular mechanisms of how these occur remains unclear. Other groups have shown that statins may trigger mevalonate-independent pathways partially via a calcium increase and p38 MAPK activation, which induces UPR and cytoprotection in RAW264.7 cells 29 . Further complicating the story is the involvement of various inflammatory molecules, such as tumor necrosis factor (TNF) and interleukin 8 (IL-8) in all three treated groups and interleukin 20  in the toxin treated groups. TNF is a pro-inflammatory cytokine leading to the activation of one if not all of three pathways: nuclear factor kappa-light-chain-enhancer of activated B cells (NF-kappaB) activation, activation of MAPK pathways and induction of apoptosis 30 . Potentially as a consequence of the activation of NF-kappaB, IL-8 production increased. IL-8 production has been previously seen in response to various bacterial toxins, such as Staphylococcal α toxin 31 , aerolysin 32 , pneumolysin 33 and Listeria monocytogenes toxin listeriolysin O 34 . Statins themselves have been shown to have both proinflammatory and anti-inflammatory properties. It has been shown that lipophilic statins (i.e. fluvastatin, simvastatin, atorvastatin, and lovastatin) were able to induce proinflammatory responses in both LPS stimulation 35 and Mycobacterium tuberculosis exposure 36 . Interleukin 20, a member of the interleukin-10 family of cytokines, is a pro-inflammatory cytokine that is associated with both psoriasis (skin inflammation) and atherosclerosis (a chronic inflammatory disease resulting in lipid deposits). Interestingly, IL-20 has been previously shown to play a rather negative role in host defense; infection with Staphylococcus aureus led to the early up regulation of IL-20 family members (IL-19, IL-20 and IL-24), inhibiting the local generation of IL-1β and IL-17A. This ultimately leads to an increase in infection severity and led the authors to hypothesize that S. aureus may express a virulence factor that induces this expression 22 . Our study is the first case showing this same up regulation of IL-20 in response to the pore-forming toxin from S. pneumonia. Note: IL-24 was also up regulated in any toxin-exposed cells but not at a level that was statistically significant. Statin co-treatment did not decrease the expression of these IL-20 family members and its protective effects must be due to a different mechanism.

KEGG Canonical Pathway
Furthermore, a role for phospholipases in statin conferred cellular protection is hinted at in this dataset. A2 phospholipases are a class of enzymes that catalyzes the release of arachidonic acid and lysophospholipids from phospholipids. Arachidonic acid is then further modified by cyclooxygenases into eicosanoids, including prostaglandins and leukotrienes, which can be both anti-inflammatory and inflammatory mediators depending on other factors 37 . A2 phospholipases are expressed and released (only the secreted subtype) by many human immune cells: macrophages, monocytes, T cells, mast cell and neutrophils. Many different isoforms of phospholipases are increased in concentration in the blood of patients with inflammatory or autoimmune diseases 38,39 . In normal systems, secreted phospholipases have been shown to be important in fighting both viral and bacterial infections, of which group IIA is the best characterized. GIIA sPLA2 has displayed antibacterial activity towards Gram-positive bacteria including S. aureus and L. monocytogenes 40,41 . Although sPLA2-IIA was not pulled out in this dataset, it has been shown by another group to be induced by statins in a dose dependent manner 42 . Overexpression of sPLA2-IIA in transgenic mice leads to improved clearance of bacteria in the lung and a decrease in mortality from Staphylococcus aureus infection 43 . This function is calcium-dependent and negated in the presence of EGTA, which is consistent with our BAPTA-AM data 44,45,46 . It is thought that the antibacterial activity depends on the ability of PLA2 to hydrolyze the bacterial cell membrane, due to the highly charged surface in the case of sPLA2-IIA 40 . Other members of the phospholipase A2 family have been shown to have antibacterial activity, with GX family members (PLA2G10 is of this family) being the second most potent sPLA2 against gram-positive bacteria (sPLA2-IIA is the most potent) 41 .
Since statins are competitive inhibitors of a key enzyme regulating cholesterol biosynthesis, it was logical to investigate a potential global lipid metabolic shift. According to the LIPID Metabolites And Pathways Strategy (LIPID MAPS), there are 8 main categories of lipids: Fatty acyls/acids (FA), Glycerolipids (GL), Glycerophospholipids (GP), Sphingolipids (SP), Sterol lipids (ST), Prenol Lipids (PR), Saccharolipids (SL) and Polyketides (PK). A majority of the lipids found to be statistically significant in our dataset were members of the glycerophospholipids category. Glycerophospholipids are key components of the lipid bilayer of cells and the following types are found in biological membranes: phosphatidylcholine (PC or lecithin), phosphatidylethanolamine (PE), phosphatidylinositol (PI) and phosphatidylserine (PS). We have hypothesized a role of statins in membrane remodeling and the fact that we have pulled out key lipids in membrane structures reinforces this idea. Interestingly, PLY exposure to the cells alone led to a dramatic shift in the saturated PC lipid profile (top statistically significant hits: PC 32:1, PC 34:2, and PC 36:1), the most abundant lipid subtype found in the plasma membrane. We found that statins appear to affect other lipid subtypes, such as triglycerides (TG), in addition to PC lipids. TG 54:4, one of the statistically significant TG lipids from the statin treated group, has been shown by another group to be up regulated in response to inflammation in a mouse macrophage LPS induced inflammatory model and expression of this lipid did not change with statin treatment 47 . Another subtype unique to the airway epithelium was the change seen in phosphatidylinositol (PI) lipids. Mouse macrophages do not show any change in expression with statin treatment whereas our work shows an increase in both PI 36:3 and PI 38:4 in response to simvastatin treatment. This is not surprising, due to the fact that PI lipids are more common in the airway than in macrophages 48 . Similarly, in the mouse macrophage, they showed that statin treatment resulted in a blockade of the formation of sterol esters, pushing temporal changes in fatty acids in response to inflammation to other categories, such as glycerolipids, glycerophospholipids and/or sphingolipids, which were all late onset lipids 47 . Further reinforcing our hypothesis for a role of phospholipases in statin-dependent protection, this group also saw an increase in expression of various PLA2 family members 47 .
The lipid complexity was also affected by simvastatin. Under statin exposure, an increase in unsaturated PC and PI lipids was detected (i.e. PC 37:4 and PC 35:4). Chain length and unsaturation levels affect both membrane thickness and membrane fluidity, both of which can affect cellular membrane integrity and ultimately, host cellular defense 49 . Looking directly at our work with lysoPC(18:1), our results showing a slight increase in cell viability with lipid co-treatment are counter to our initial expectations that the reduced lysoPC(18:1) conferred the cellular protection against PLY. However, since lysoPC(18:1) is known to increase during the initial stages of apoptosis and this lipid was significantly decreased when simvastatin was given before toxin exposure, our results may indicate that the decrease of lysoPC(18:1) could be a downstream effect caused from reduced cell death triggered by PLY intoxication. This is due to the ability of lysoPC lipids to interact with various pro-apoptotic proteins of the bcl-2 family, having been shown to decrease the amount of bid protein associated with mitochondria [50][51][52] . It will be interesting to see what role lipid dynamics play in cellular protection as this field advances with the advent of more commercially available lipids.
A possible explanation for the lipid metabolic and complexity shift under simvastatin exposure could be due to simvastatin-regulated genes. Comparing our dataset with another RNA-seq dataset involving the exposure of human liver cells (HepG2 cell line) to atorvastatin, using the same log2FC cutoff of 1.2 used by this group 53 , we found 22 transcripts in common between these two datasets. Many of the genes enriched, such as mvd and hmgcs1, are involved in mevalonate and cholesterol biosynthesis pathways. It is plausible that statins in general up-regulate mevalonate pathway genes to compensate for a reduced output from the mevalonate pathway in both cell types. On the other hand, the rest of the regulated genes could be potential unique statin-mediated response in either cell type.
Overall, we have shown a role for statins in protection against the negative effects of pore-forming toxins. Through the use of next-generation sequencing and lipidomics, we have provided global insights into how this protection may occur. Furthermore, our findings reinforce the idea that statin use may be potentially developed as an adjuvant therapy to pore-forming toxin-related bacterial infection in the lung.

Methods
Ethics Statement. Excised human airways were collected from the NIH-sponsored National Disease Research Interchange (NDRI, Philadelphia, PA) procurement through informed consent of the donors or next-of-kin. The investigators performing the experiments are not involved in the recruitment or direct contact with the donors and their relatives. In addition, the tissues were anonymized and the investigators have no access to the personal information of the donors. The use of cells isolated from those airways has been approved by Institute Human Subject Review Committee of University of California Davis. The methods were carried out in accordance with approved guidelines.
Cell Culture and Reagents. Simvastatin was purchased from Sigma-Aldrich. C18:1 lysophosphatidylcholine was purchased from Cayman. ROCK inhibitor Y-27632 was purchased from Reagents Direct. 6xHis tag-fused pneumolysin was expressed in and purified from an Escherichia coli M15 (pREP4) 54 . Residual LPS was removed by passage over AffinityPak Detoxi-Gel (Endotoxin Removal Gel; Thermo Fisher Scientific) following the manufacturer's instructions.
Illumina library preparation/sequencing. NHBE cells were seeded at 5 × 10 5 cells per well of 6-well plates and grown overnight. Three wells were used per treatment. In two days, cells were exposed to simvastatin (1 μ M, Sigma-Aldrich) or vehicle control in BEGM medium for 24 hours. The exposed cells were then treated with pneumolysin (400 ng/mL) or mock control for 4 hours. This resulted in 4 different classes from three different patients: vehicle control (N), only statin (S), only toxin (T), and both statin and toxin (B). Total RNA was extracted and purified using RNeasy kit (Qiagen) per manufacturer's instructions. The RNA integrity (RIN Score > 9.0) and quantity was determined on the Agilent 2100 Bioanalyzer (Agilent) using the RNA 6000 nano kit following manufacturer's instructions. RNA samples were prepared by Beijing Genomics Institute using the Illumina TruSeq RNA sample preparation kit (Illumina). Briefly, total RNA (2 μ g) was used for poly-A mRNA selection using oligodT magnetic beads, with two rounds of enrichment. The samples were then fragmented using both heat and high salt conditions. cDNA synthesis, multiplex indexing, and PCR amplification was then completed accordingly, resulting in a dsDNA library for each of the 12 samples. Indexed libraries were pooled in equal concentrations and run on a HiSeq2000 for SR50. HiSeq2000 reads were aligned to the Hg19 transcriptome using Burrows-Wheeler Transform allowing for 4% mismatch followed by conversion to SAM files using SAMtools 55 . Indexed samples were then merged at the SAM file step, resulting in 12 samples. The read counts were extracted with a perl script (provided upon request), and then imported into R. Data analysis was completed using edgeR 56 and network figures were created using cytoscape 57 with the iRegulon plugin 25 . The raw data files have been uploaded to NCBI, under the following bioproject ID: PRJNA277544.
Gene Set Enrichment Analysis (GSEA). GSEA was run using the graphical user interface javaGSEA desktop application, downloaded from the GSEA website (http://www.broadinstitute.org/gsea). Aligned data was adjusted to HUGO gene name format in order to upload into GSEA, using raw reads as expression values. A class file was created to indicate the group of each of the 12 samples to the GSEA program.
Lipidomics. HBE1 cells were seeded at 5 × 10 5 cells per well of 6-well plates and grown overnight.
Three wells were used per replicate of each treatment. In two days, cells were exposed to simvastatin (1 μ M) or vehicle control in growth medium for 24 hours. After medium change, the pre-exposed cells were then treated with pneumolysin (400 ng/mL) or mock control for four hours. This resulted in 4 different classes with seven replicates per condition: vehicle control (EC), only statin (SC), only toxin (EP), and both statin and toxin (SP). Cells were washed three times with PBS and pelleted at 10 RCF for 10 minutes at 4 °C. The PBS was removed and the pellet was snap frozen and stored at − 80 °C until used. Samples were thawed and extracted with MTBE as described previously 58 . Briefly, the vortexed/grinded homogenized sample (30 μ L) was placed in a 1.5 mL vial. Cold methanol (225 μ L) and MTBE (750 μ L) were added to each sample. Samples were vortexed and shaken for six minutes at 4 °C. Distilled water was added to bring total volume to 1.2 mL and samples were centrifuged for 2 minutes at 14,000 RCF. Samples were dried to complete dryness and reconstituted in with MeOH/chloroform (100 μ L, 9/1, v/v). For mass spectrometric analysis, sample (10 μ L) was then added to MeOH/chloroform (90 uL) containing 7.5 mM ammonium acetate in a well of a 96 well plate and then sealed with aluminum foil. Mass spectrometric analysis was performed on a TSQ Quantum Ultra Plus triple-quadrupole mass spectrometer (Thermo Fisher Scientific) equipped with an automated nanospray apparatus (Nanomate HD, Advion Bioscience Ltd.). MS survey scans were acquired in both positive and negative ion mode. These scans were processed with GeneData Expressionist Refiner MS software (GeneData) and further validated using the NIST MS Search program to compare MS/MS data with an in-house library called LipidBlast 59 . Intensities were further analyzed using Statistica 9.0 software (StatSoft). Heatmap was created with the assistance of R 3.0.1 (http://www.r-project.org/), using Euclidean distance measures and the Ward clustering algorithm.
Cytotoxicity assay. HBE1 cells were seeded at 10,000 cells per well of 96-well plates and grown overnight. Cells were exposed to simvastatin, C18:1 lysophosphatidylcholine at indicated concentrations or vehicle control in growth medium for 24 hours. The exposed cells were then treated with pneumolysin at indicated concentrations for 4 hours and cell viability was analyzed with CellTiter-Glo® Luminescent Cell Viability Assay (Promega) following recommended protocol.
Real-time RT-PCR. The procedures for RNA extraction and cDNA generation were described in our previous study 60 . Maxima SYBR Green qPCR Master Mix no Rox (Thermo Scientific) was used to perform Real-time PCR reactions on an Applied Biosystems 7900HT. Expressions were normalized to GAPDH for all reactions. Sequences of primers used in this study are included in Table ST3.