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.

Due to the constant exposure to the outside world, the lung presents a prime target for bacterial and viral infections. For the normal healthy adult airway, mucociliary clearance and acute immune responses usually work well enough to decrease the incidence of infection. Yet, in the compromised airway, such as seen in asthmatics, this is insufficient1. One of the most frequent bacterial infections found in asthmatics is Streptococcus pneumoniae, or pneumococcus2. This bacterial infection is also found in the general public and is the most common cause of bacterial pneumonia, as well as the root cause of some cases of bacterial meningitis, bacteremia, and otitis media3. S. pneumoniae infection is associated with both high morbidity and high mortality rates34. Although exact statistics are hard to determine, the World Health Organization estimates that 1.6 million deaths were caused by S. pneumoniae infection in 2008 (world-wide), with 333,000 to 529,000 of these occurring in children younger than 5 years56.

Pneumococcus is a polysaccharide capsulated gram-positive bacteria and more than 90 serotypes have been identified that vary in virulence and drug resistance34. Like many other bacteria, S. pneumoniae 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 pathogens7. Pneumolysin is secreted as a soluble monomer, which upon binding to host cell membranes via cholesterol, can oligomerize to form β-barrel pores and enter cells89. 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 viability10. 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 biosynthesis11. 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 mortality12. 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 users13. In contrast, laboratory studies using animal models of infection have shown protection against bacterial infections under statin use, mainly focusing on Staphylococcus aureus141516. 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 composition11. Currently, there are only a limited amount of reports involving the lipidomic analysis of statins, with a majority carried out with human plasma samples171819. 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. lysoPC 18:1, PE 12:0) while decreasing longer and unsaturated lipids (i.e. PC 36:2, PG 36:1, PS 38:2, PS 39:8, PS 40:2). In contrast, statin pretreatment led to a shift towards unsaturated lipids, as seen in PI 36:3, PI 38:4, and longer chained lipids (i.e. TG 52:4, TG 54:4). This lipid data is shown in figs 2 and 3.

Figure 1

Global identification of lipid members in human bronchial cells (HBE1). Heatmap showing hierarchical clustering of 89 identified lipid members via mass spectrometry in the 4 classes of treatments: 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). Each experimental exposure was completed seven times in HBE1 cells and all independent samples are shown as z-score data. Intensities were normalized by class and underwent z scale centering by individual lipid. Clustering analysis was completed using Euclidean distance measures and the Ward algorithm.

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).

Lipid	p Value	Control vs Toxin	Control vs Statin	Control vs Both	Toxin vs Both	Statin vs Both	Statin vs Toxin
plasmenyl-PC 38:4	0.000002	ns	***	***	*	ns	***
PC 30:0	0.000004	***	***	*	ns	**	ns
PC 32:1	0.000006	***	ns	ns	***	*	ns
PA 34:2	0.00001	***	***	ns	ns	**	ns
PC 32:0	0.000016	***	**	ns	***	**	ns
PC 34:2	0.000018	***	***	ns	**	**	ns
PC 36:1	0.000038	***	*	ns	***	*	ns
PC 36:2	0.000078	***	*	ns	***	*	ns
PS 40:2	0.000088	*	ns	ns	*	ns	***
PS 38:2	0.000093	ns	*	ns	**	ns	***
TG 52:4	0.000242	ns	**	**	*	**	*
PC 37:4	0.000566	ns	**	**	ns	ns	*
PC 35:4	0.000694	ns	***	*	ns	ns	*
lysoPC 18:1	0.001047	ns	ns	ns	**	ns	**
plasmenyl-PC 38:5	0.001055	**	**	ns	ns	ns	ns
TG 52:3	0.001548	ns	ns	ns	**	**	**
PC 36:2	0.001665	***	*	ns	***	*	ns
TG 54:4	0.001816	ns	**	**	ns	ns	ns
PG 36:1	0.003718	***	*	ns	***	*	ns
PS 36:1	0.004565	ns	ns	ns	ns	ns	*
PS 39:8	0.004883	*	ns	ns	ns	ns	*
PC 30:1	0.00537	ns	ns	**	*	ns	ns
PS 32:1	0.005472	ns	ns	**	ns	ns	ns
PG 34:1	0.007567	ns	ns	ns	ns	ns	*

Figure 2

Pore-forming toxin exposure led to a decrease in lipid complexity, resulting in more saturated, short chain lipids. Six box plots characteristic of this trend are shown for the following lipids: PC 36:2, PE 12:0, PG 36:1, PS 38:2, PS 39:8 and PS 40:2. 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.

Figure 3

Statin pretreatment led to an increase in cell survival and lipid complexity, resulting in more unsaturated, long chain lipids. A.) Five box plots characteristic of this trend are shown for the following lipids: PI 36:3, PI 38:4, TG 52:4, TG 54:4 and lysoPC 18:1. 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. 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.

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 enzymes20. 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 pecies">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 2

Top twenty differentially expressed genes from the mRNA-seq data, as compared to no treated controls. mRNA-seq reads were normalized using scaling factors for library size. P values were determined via a generalized linear model (GLM) likelihood ratio test within the edgeR package.

Gene Symbol	Gene Definition	Log2 Fold Change, Treatment/None	p Value
TOXIN
EGR1	Early growth response 1	5.0821124356	1.02E-05
EGR2	Early growth response 2	4.5849534080	5.40E-06
IL20	Interleukin 20	4.5784428610	0.0002824629
C11orf96	Chromosome 11 open reading frame 96	4.4944040287	4.75E-06
ATF3	Activating transcription factor 3	4.4694850002	8.68E-06
TNF	Tumor necrosis factor	3.8800929090	0.0004175837
GEM	GTP binding protein overexpressed in skeletal muscle	3.4577584126	3.08E-05
NR4A1	Nuclear receptor subfamily 4, group A, member 1	3.2596075994	8.62E-07
CXCL2	Chemokine (C-X-C motif) ligand 2	3.2279840457	8.24E-05
IL8	Interleukin 8	3.2133483633	7.24E-06
NR4A3	Nuclear receptor subfamily 4, group A, member 3	3.0059522416	0.0078231198
FOSB	FBJ murine osteosarcoma viral oncogene homolog B	2.8638213316	0.0002123400
CREBRF	CREB3 regulatory factor	2.2232930679	9.67E-05
GADD45B	Growth arrest and DNA-damage-inducible, beta	2.2001574628	1.39E-08
GDF15	Growth differentiation factor 15	2.1755167968	6.21E-05
RND1	Rho family GTPase 1	2.0872777098	1.49E-08
PPP1R15A	Protein phosphatase 1, regulatory subunit 15A	2.0412893595	2.04E-05
BIRC3	Baculoviral IAP repeat containing 3	1.9263031522	0.0018212122
MAP2K6	Mitogen-activated protein kinase kinase 6	1.8155637257	0.0086922392
FOS	FBJ murine osteosarcoma viral oncogene homolog	1.7937408311	2.40E-07
STATIN
FXYD3	FXYD domain containing ion transport regulator 3	4.433078947	1.45E-05
C11orf96	Chromosome 11 open reading frame 96	3.72461041	0.000145175
TNF	Tumor necrosis factor	3.695235714	0.001722366
NFASC	Neurofascin	3.635372308	8.42E-05
C6orf223	Chromosome 6 open reading frame 223	3.528111821	0.001664323
CD9	CD9 molecule	3.305106179	0.007018955
MAFF	V-maf avian musculoaponeurotic fibrosarcoma oncogene homolog F	3.147336747	3.14E-06
ATF3	Activating transcription factor 3	3.126219473	0.002610577
CXCL2	Chemokine (C-X-C motif) ligand 2	3.034548443	6.97E-06
SLC27A2	Solute carrier family 27 (fatty acid transporter), member 2	2.898860433	0.001683999
HSD17B7	Hydroxysteroid (17-beta) dehydrogenase 7	2.824996536	9.03E-06
REC8	REC8 meiotic recombination protein	2.793153813	1.23E-05
SLC29A2	Solute carrier family 29 (equilibrative nucleoside transporter), member 2	2.688913638	0.003246542
IL8	Interleukin 8	2.4637838	0.000248558
TM7SF2	Transmembrane 7 superfamily member 2, Delta(14)-sterol reductase	2.41216306	0.000137761
LGALS7B	Lectin, galactoside-binding, soluble, 7B	2.338793329	0.000115045
FDPS	Farnesyl diphosphate synthase	2.245113522	0.005511471
ACACB	Acetyl-CoA carboxylase beta	2.18279627	0.007533777
GPX2	Glutathione peroxidase 2	2.143697892	0.008187339
SSBP2	Single-stranded DNA binding protein 2	2.136914726	1.60E-05
STATIN+TOXIN
EGR2	Early growth response 2	9.003110027	8.48E-23
C11orf96	Chromosome 11 open reading frame 96	8.51357552	1.45E-14
C6orf223	Chromosome 6 open reading frame 223	7.494553998	1.30E-10
ATF3	Activating transcription factor 3	7.208837244	8.96E-11
EGR1	Early growth response 1	7.163627625	7.66E-09
NR4A1	Nuclear receptor subfamily 4, group A, member 1	6.922380993	7.32E-20
IL8	Interleukin 8	6.778865409	5.51E-16
IL20	Interleukin 20	6.767729321	1.19E-07
TNF	Tumor necrosis factor	6.696143663	2.67E-08
FOSB	FBJ murine osteosarcoma viral oncogene homolog B	5.910606615	9.77E-12
CXCL2	Chemokine (C-X-C motif) ligand 2	5.589639741	3.35E-10
GEM	GTP binding protein overexpressed in skeletal muscle	5.479633585	4.96E-10
SOCS3	Suppressor of cytokine signaling 3	4.74582123	1.18E-17
ZNF460	Zinc finger protein 460	4.462020957	8.90E-12
PLA2G3	Phospholipase A2, group III	4.35384952	3.26E-06
TNFAIP3	Tumor necrosis factor, alpha-induced protein 3	4.253405429	3.07E-17
BIRC3	Baculoviral IAP repeat containing , apoptosis inhibitor	4.24819277	6.31E-10
KLF2	Kruppel-like factor 2	4.160537722	1.09E-08
RASGEF1B	RasGEF domain family, member 1B	4.15890662	8.45E-17
CCL20	Chemokine (C-C motif) ligand 20	4.067584886	1.09E-09

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 protein response (UPR), an endoplasmic reticulum stress response21. 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 pathways22. 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 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 viruses23.

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 synthesis24. 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 nucleus24. 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 iRegulon25 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.

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.

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 Enrichpecies">ment 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 dataset26. 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.

Table 3

Enriched Gene Ontology (GO) pathways in the transcriptome of treatment groups as compared to no treated controls. Nominal p values are calculated with the GSEA program by using an empirical phenotype-based permutation test. False discovery rate is adjusted for using an adapted Benjamini–Hochberg procedure.

GO Canonical Pathway	Number of Genes in Set	NOM p-val	FDR q-val
TOXIN
RESPONSE TO EXTRACELLULAR STIMULUS	22	0	0.002938457
NITROGEN COMPOUND CATABOLIC PROCESS	18	0	0.007242801
CATION TRANSPORT	61	0	0.019453308
RESPONSE TO NUTRIENT LEVELS	18	0	0.019539576
INFLAMMATORY RESPONSE	65	0	0.020307241
TRANSLATION	151	0	0.024203578
STEROID METABOLIC PROCESS	43	0	0.02606539
MONOCARBOXYLIC ACID METABOLIC PROCESS	65	0	0.039081186
GLYCEROPHOSPHOLIPID BIOSYNTHETIC PROCESS	29	0	0.04601446
LIPID METABOLIC PROCESS	227	0	0.04833658
STATIN
MONOCARBOXYLIC ACID METABOLIC PROCESS	65	0	0
LIPID METABOLIC PROCESS	227	0	0.0018675
GLUCOSE METABOLIC PROCESS	21	0	0.002207291
FATTY ACID METABOLIC PROCESS	50	0	0.002241
CELLULAR LIPID METABOLIC PROCESS	185	0	0.002287963
STEROID METABOLIC PROCESS	43	0	0.003534764
FATTY ACID OXIDATION	16	0	0.006060679
CARBOXYLIC ACID METABOLIC PROCESS	135	0	0.006493585
LIPID BIOSYNTHETIC PROCESS	80	0	0.006537986
INFLAMMATORY RESPONSE	65	0	0.03650236
STATIN+TOXIN
LIPID METABOLIC PROCESS	227	0	0.007665763
CATION TRANSPORT	61	0	0.021874871
GLUCOSE METABOLIC PROCESS	21	0	0.023219533
CELLULAR LIPID METABOLIC PROCESS	185	0	0.030281333
RESPONSE TO EXTERNAL STIMULUS	164	0	0.032070216
MONOCARBOXYLIC ACID METABOLIC PROCESS	65	0	0.03301661
RESPONSE TO NUTRIENT LEVELS	18	0.008547009	0.04540056
DEFENSE RESPONSE	126	0	0.04604433
MEMBRANE LIPID BIOSYNTHETIC PROCESS	44	0.007407407	0.047668364
INFLAMMATORY RESPONSE	65	0	0.048929326

Table 4

Enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in the transcriptome of treatment groups as compared to no treated controls. Nominal p values are calculated with the GSEA program by using an empirical phenotype-based permutation test. False discovery rate is adjusted for using an adapted Benjamini–Hochberg procedure.

KEGG Canonical Pathway	Number of Genes in Set	NOM p-val	FDR q-val
TOXIN
RIBOSOME	85	0	0
GLUTATHIONE METABOLISM	41	0	0
LYSOSOME	110	0	0
BIOSYNTHESIS OF UNSATURATED FATTY ACIDS	18	0	0.008914938
STEROID BIOSYNTHESIS	15	0	0.013855341
LEISHMANIA INFECTION	43	0	0.014490775
ARACHIDONIC ACID METABOLISM	34	0	0.015463158
CARDIAC MUSCLE CONTRACTION	38	0	0.015587258
STEROID HORMONE BIOSYNTHESIS	30	0	0.017814009
METABOLISM OF XENOBIOTICS BY CYTOCHROME P450	42	0.00990099	0.021507936
STATIN
STEROID BIOSYNTHESIS	15	0	0
BIOSYNTHESIS OF UNSATURATED FATTY ACIDS	18	0	0
ECM RECEPTOR INTERACTION	50	0	0
METABOLISM OF XENOBIOTICS BY CYTOCHROME P450	42	0	0
ADIPOCYTOKINE SIGNALING PATHWAY	52	0	0.001815
GLUTATHIONE METABOLISM	41	0	0.002716081
DRUG METABOLISM CYTOCHROME P450	39	0	0.005934439
STEROID HORMONE BIOSYNTHESIS	30	0	0.006539325
DILATED CARDIOMYOPATHY	44	0	0.007513437
PYRUVATE METABOLISM	31	0	0.008348262
STATIN+TOXIN
RIBOSOME	85	0	0
STEROID BIOSYNTHESIS	15	0	0
BIOSYNTHESIS OF UNSATURATED FATTY ACIDS	18	0	0
METABOLISM OF XENOBIOTICS BY CYTOCHROME P450	42	0	0
GLUTATHIONE METABOLISM	41	0	0
STEROID HORMONE BIOSYNTHESIS	30	0	0
LYSOSOME	110	0	0.000965857
ADIPOCYTOKINE SIGNALING PATHWAY	52	0	0.009771074
PPAR SIGNALING PATHWAY	40	0	0.011506875
ARACHIDONIC ACID METABOLISM	34	0	0.011509539

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 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. pneumoniae27. 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 production28. 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 cells29. 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 (IL-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 apoptosis30. 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 α toxin31, aerolysin32, pneumolysin33 and Listeria monocytogenes toxin listeriolysin O34. 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 stimulation35 and Mycobacterium tuberculosis exposure36. 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 expression22. 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.

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 factors37. 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 diseases3839. 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. monocytogenes4041. 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 manner42. Overexpression of sPLA2-IIA in transgenic mice leads to improved clearance of bacteria in the lung and a decrease in mortality from Staphylococcus aureus infection43. This function is calcium-dependent and negated in the presence of EGTA, which is consistent with our BAPTA-AM data444546. 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-IIA40. 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 treatment47. 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 macrophages48. 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 lipids47. 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 members47.

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 defense49. 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 mitochondria505152. 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 group53, 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 pecies">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.

HBE1 cells, a human papillomavirus immortalized bronchial epithelial cell line and a gift from JR Yankaskas at University of North Carolina, were cultured using BEGM Bronchial Epithelial Cell Growth Medium (Lonza, Basel, Switzerland), minus the retinoic acid: Ham’s F12/Dulbeco’s Modified Eagle (DME) medium (1:1) supplemented with insulin (5 μg/mL), transferrin (5 μg/mL), epidermal growth factor (10 ng/mL), dexamethasone (0.1 μM), cholera toxin (10 ng/mL), bovine hypothalamus extract (15 μg/mL), and bovine serum albumin (0.5 mg/mL). Regarding Normal Human Bronchial Epithelial (NHBE) cells, protease-dissociated tracheal and primary bronchial epithelial cells were collected from tissues of patients without diagnosed lung-related diseases and plated in BEGM Bronchial Epithelial Cell Growth Medium (Lonza) minus the retinoic acid with the addition of 5 μM ROCK inhibitor. All cells were maintained at 37 °C in a humidified incubator with 5% CO2.

Illumina library preparation/sequencing

NHBE cells were seeded at 5 × 105 cells per well of 6-well plates and grown overnight. Three wells were used per treatpecies">ment. 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 SAMtools55. 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 edgeR56 and network figures were created using cytoscape57 with the iRegulon plugin25. 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 × 105 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 previously58. 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 LipidBlast59. 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 study60. 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.

Additional Information

How to cite this article: Statt, S. et al. pan class="Chemical">Lipidome and Transcriptome Profiling of Pneumolysin Intoxication Identifies Networks Involved in Statin-Conferred Protection of Airway Epithelial Cells. Sci. Rep. 5, 10624; doi: 10.1038/srep10624 (2015).