Campylobacter spp. are one of the most important food-borne pathogens, which are quite susceptible to environmental or technological stressors compared to other zoonotic bacteria. This might be due to the lack of many stress response mechanisms described in other bacteria. Nevertheless, Campylobacter is able to survive in the environment and food products. Although some aspects of the heat stress response in Campylobacter jejuni are already known, information about the stress response in other Campylobacter species are still scarce. In this study, the stress response of Campylobacter coli and Campylobacter lari to elevated temperatures (46°C) was investigated by survival assays and whole transcriptome analysis. None of the strains survived at 46°C for more than 8 h and approximately 20% of the genes of C. coli RM2228 and C. lari RM2100 were differentially expressed. The transcriptomic profiles showed enhanced gene expression of several chaperones like dnaK, groES, groEL, and clpB in both strains, indicating a general involvement in the heat stress response within the Campylobacter species. However, the pronounced differences in the expression pattern between C. coli and C. lari suggest that stress response mechanisms described for one Campylobacter species might be not necessarily transferable to other Campylobacter species.

Introduction

Campylobacter is one of the most common causative agents of bacterial food-borne gastroenteritis in humans worldwide. The campylobacteriosis (with clinical symptoms like diarrhea, abdominal cramps, and fever) is mostly self-limiting but sequelae like Guillian–Barré syndrome and reactive arthritis have been described (Poropatich et al., 2010). Handling and consumption of raw or undercooked meat has been identified as a main source for human infections (Alter et al., 2011). Campylobacter belongs to the heterogeneous class of Epsilonproteobacteria and the species mostly detected in diseased humans are Campylobacter jejuni followed by Campylobacter coli and to a lesser extent Campylobacter lari and Campylobacter upsaliensis (Whiley et al., 2013). Whereas C. jejuni shows high prevalences in poultry and cattle, C. coli is usually associated with pigs and C. lari is frequently detected in shellfish and shorebirds (Miller et al., 2008; EFSA and ECDC, 2014). However, in these animals, Campylobacter belongs to the commensal microbiota. By comprehensive genome comparisons of several strains, it has been shown that C. jejuni is more closely related to C. coli compared to C. lari (Zhou et al., 2013).

Campylobacter spp. possess fastidious growth requirements and are less resistant against environmental and technological stressors compared to other zoonotic pathogens (Park, 2002). This might be explained by the absence of several stress response regulators typically involved in the regulation of stress response to various stressors in other Gram-negative bacteria like Escherichia coli or Salmonella, e.g., alternative sigma factors like RpoS (σ³⁸) as well as other transcription regulators such as CspA, Lrp, SoxRS, and OxyR (Alter and Scherer, 2006). Unlike E. coli, which harbors seven sigma factors, only RpoD (σ⁷⁰), RpoN (σ⁵⁴), and RpoF/FliA (σ²⁸) were found in Campylobacter spp. (Fouts et al., 2005). However, Campylobacter spp. is able to survive in the environment and to overcome the barriers along the food chain. This suggests the existence of alternative regulatory mechanisms or a wider role of known regulatory factors. For example, the transcription of both the peroxide as well as the superoxide defense genes are regulated by PerR in C. jejuni while in E. coli and Salmonella OxyR regulates the expression of the peroxide defense regulon and SoxR of the superoxid defense regulon (Kim et al., 2015).

Several authors have described the response of C. jejuni to increased temperatures. Konkel et al. (1998) were able to identify 24 proteins preferentially synthesized following heat stress at 46°C. One of these proteins was identified as the DnaJ chaperon, shown to be necessary for survival at elevated temperatures. Further proteome analyses showed 18 differentially expressed proteins induced by increasing temperature from 37 to 42°C (Zhang et al., 2009). The gene expression profiles associated with these responses have also been investigated by microarray analysis (Stintzi, 2003). This study revealed an increased expression of common genes encoding heat stress proteins like chaperones and proteases (e.g., groELS, grpE, dnaK, dnaJ, clpB, lon) and membrane associated proteins (e.g., galE, gmhA2). In addition, the chaperon activity of the serine protease HtrA is involved in the heat stress response of C. jejuni (Baek et al., 2011). The alternative sigma factor σ³² (RpoH) mediating the expression of heat stress related genes in E. coli is missing in Campylobacter (Parkhill et al., 2000; Yura et al., 2000). So far, RacR, HspR, and HrcA have been identified as regulators of the heat stress response in C. jejuni (Holmes et al., 2010; Apel et al., 2012), while RpoN (σ⁵⁴), known to be involved in various stress responses, is not involved in the heat stress response of C. jejuni (Hwang et al., 2011).

Survival strategies of the related species C. coli and C. lari at high temperatures are largely unexplored. Thus, our study aimed to compare the effects of heat stress (46°C) on the two Campylobacter species C. coli and C. lari by (i) survival studies and (ii) whole transcriptome (RNA-seq) analyses, to identify genes involved in the heat stress response.

Materials and Methods

Bacterial Strains, Media, and Growth Conditions

Campylobacter strains were grown on Mueller–Hinton agar containing 5% sheep blood (MHB; OXOID, Wesel, Germany) for 48 h or Brucella broth for 24 h (BB; BD, Heidelberg, Germany) at 37°C in microaerobic conditions (6% O₂, 7% CO₂, 7% H₂, 80% N₂) generated by the Mart Anoxomat system (Drachten, Netherlands). The whole genome sequenced strains C. coli RM2228 and C. lari RM2100 were selected according to their clinical importance. The C. coli strain RM2228, isolated from a chicken carcass, belongs to the phylogenetic clade1, which are mainly responsible for human infections (Sheppard et al., 2010). The C. lari strain RM2100 was isolated from a child with watery diarrhea (Miller et al., 2008). Two further field strains of each species were included for survival assays. All strains are listed in Supplementary Table S5.

Survival Studies

For survival studies, pre-cultures were diluted in BB to an optical density of 0.01 at 600 nm (approximately 7 log₁₀ CFU/ml) and incubated at 46°C in microaerobic conditions in a static cultivation mode. The cell numbers were determined over a 24 h period by plating serial dilutions on MHB. These plates were incubated at 37°C for 48 h and cell counts shown as log₁₀ CFU/ml. Experiments were performed at three individual time points with technical duplicates.

Quantitative Real-Time PCR

The expression of the selected genes clpB, grpE, dnaK, groEL, groES, cbpA, and dnaJ of C. coli and C. lari were analyzed over 60 min stress exposure, induced by increasing temperatures from 37°C to 46°C. RNA extractions of three individual cultures were included for both species. Total RNA was isolated from approximately 9 log CFU using a peqGOLD Bacterial RNA Kit (Peqlab, Erlangen, Germany). Removal of genomic DNA was performed in a total volume of 40 μl containing 4 U DNase I, 40 U Ribolock, 1x DNase buffer (all Fermentas, Leon-Rot, Germany), and 28 μl of RNA. After an incubation for 15 min at 37°C, DNase was inactivated by adding 4 μl 50 mM EDTA and heating at 65°C for 10 min. First-strand cDNA was synthesized of 1 μg RNA using the RevertAid Premium First Strand cDNA synthesis kit and random hexamer primers (all Fermentas) according to manufacturer’s instructions.

Primers were designed using the Primer3 web interface¹ based on the C. coli RM2228 and C. lari RM2100 genome sequences (with primer length 17–25 nt, amplicon size 60–150 nt, and primer Tm 50–62°C). All runs were performed in a 15 μl PCR mixture containing 1 μl of a 1:10 dilution of cDNA, 0.05–0.9 μM of each primer (Supplementary Table S1) and a twofold SsoFast EvaGreen Supermix (Bio-Rad, Munich, Germany). The amplification was performed by pre-heating for 30 s at 94°C and 40 cycles of 94°C for 5 s followed by annealing for 10 s (annealing temperatures are shown in Supplementary Table S1). Quantitative real-time PCR data were processed using the CFX Manager Software (Bio-Rad) with thiC and rpoA as housekeeping genes. ΔΔCq values of all samples were determined based on two technical replicates. Specificity of the amplification product was confirmed by melting curve analysis.

RNA Isolation for Whole Transcriptome Analysis

As quantitative real-time PCR analysis indicated still enhanced mRNA-level for several genes after 30 min heat stress, this time point was chosen for whole transcriptome analysis. RNA was isolated from two independent cultures for each species. Total RNA was extracted from approximately 9 log CFU of the two species C. coli and C. lari after 30 min cultivation at 37°C and of heat stress at 46°C using a hot phenol/lysozyme method (Blomberg et al., 1990). Therefore, the bacterial culture was mixed 1:5 with stop-mix (95% v/v ethanol, 5% v/v water saturated phenol, −20°C). Cells were pelleted by centrifugation at 1500 × g, supernatant discarded, and samples stored at -80°C. Pellets were resuspended in 600 μl 0.5 mg/ml lysozyme in TE (pH 8.0), 60 μl 10% w/v SDS were added and incubated at 64°C for 2 min. Afterward 66 μl 1 M NaOAc (pH 5.2) and 750 μl Roti-Aqua phenol (Carl Roth, Karlsruhe, Germany) were added. During incubation at 64°C for 6 min, samples were mixed by inversion 6–10 times. Tubes were placed on ice and centrifuged for 15 min at 12,000 × g (4°C). Aqueous layer was transferred in a 2 ml Phase Lock Gel-Heavy tube (VWR, Dresden, Germany), 750 μl chloroform (Carl Roth) added and mixed by inversion. After centrifugation for 12 min at 12,000 × g, the aqueous layer was transferred, 1.4 ml 30:1 EtOh:3 M NaOAc (pH 6.5) added and incubated overnight at −20°C. Samples were centrifuged at 12,000 × g (at 4°C for 30 min), supernatant removed, and pellet washed with 900 μl 75% v/v ethanol (-20°C). Samples were centrifuged as mentioned above, ethanol removed, and the pellet air-dried. The RNA was dissolved by adding 100 μl nuclease-free water (Roth) and shaken for 5 min at 65°C (800–1000 r/min). DNA was removed according to the above-described method. RNA integrity and quantity were determined on an Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany) before cDNA synthesis.

Whole Transcriptome Sequencing

Vertis Biotechnologie AG (Munich, Germany) generated libraries of two independent RNA samples for each condition. The RNA samples were poly(A)-tailed by using poly(A) polymerase. The 5’PPP were removed using tobacco acid pyrophosphatase (TAP) followed by the ligation of the RNA adapter to the 5’-monophosphate of the RNA. First-strand cDNA synthesis was performed with an oligo(dT)-adapter primer and the M-MLV reverse transcriptase. The resulting cDNA was PCR-amplified to reach a concentration of 20–30 ng/μl using a high fidelity DNA polymerase. The cDNA was purified using the Agencourt AMPure XP kit (Beckman Coulter Genomics, Essex, United Kingdom) and was analyzed by capillary electrophoresis. The primers used for PCR amplification were designed for TruSeq sequencing according to the instructions of Illumina (San Diego, CA, United States). The following adapter sequences flank the cDNA inserts: TrueSeq_Sense_primer 5’-AATGATACGGCGACCACCGAGATCTACACTCTTTCCC TACACGACGCTCTTCCGATCT-3’ and TrueSeq_Antisense _NNNNNN_primer (NNNNNN = Barcode) 5’-CAAGCA GAAGACGGCATACGAGATNNNNNNGTGACTGG AGTTC AGACGTGTGCTCTTCCGATC(dT25)-3’. The combined length of the flanking sequences is 146 bases. The libraries were sequenced with an Illumina HiSeq machine with 100 cycles in single end mode.

Bioinformatical Analysis

Read Data Analysis

The resulting sequence reads were demultiplexed and the adapter sequences were removed. After that the reads in Fastq format were quality trimmed using fastq_quality_trimmer (from the FastX suite version 0.0.13²) with a quality cut-off score of 20 and converted to fasta format using fastq_to_fasta (also from the FastX suite). The read processing [poly(A) removal, size filtering (min 12 nt length), statistics generation, coverage calculation, and normalization] was performed with the RNA-analysis pipeline READemption version 0.3.4 (Forstner et al., 2014; building upon the following libraries Biopython 1.65, pysam 0.8.1, matplotlib 1.4.3, pandas 0.16.0) which used segemehl version 0.1.7 (Hoffmann et al., 2009) for the read alignment. For C. coli RM2228, the sequences with the accession number AAFL01000001.1-38.1, DQ518170, DQ518171, as well as DQ518172.1 and for C. lari RM2100 the sequences with the accession number NC_012039 as well as NC_01240.1 were used as references (all mapping statistics are shown in Supplementary Table S2). Feature wise gene-quantification was performed (also with READemption) and used for differential gene expression analysis. According to the study of Stintzi (2003) genes with log₂ fold changes of expression below -0.58 or above 0.58 and an adjusted (Benjamini–Hochberg corrected) p-value < 0.05 calculated by DESeq2 1.6.3 were defined as regulated (Love et al., 2014).

The RNA-Seq data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus (Edgar et al., 2002) and are accessible through GEO Series accession number GSE67486³.

Orthologous Mapping

Orthologs were defined, based on the nucleotide sequences, by bidirectional best-BLAST-hit search with max. e-value of 1e-6, word size of 20, and a minimal length of 60% of both query and subject nucleotide sequence. Further some genes were determined as orthologous by similarity of gene name.

Functional Grouping

Functional groups of C. coli, C. lari, and C. jejuni CDS were based on functional categories received by eggNOG4.0 database⁴ or by functional category from orthologs. CDS belonging to two or more different categories were listed in each category. Statistical significance of enrichments of functional groups was calculated with Fisher’s exact test using GraphPad Prism version 6.07 for Windows (GraphPad Software, La Jolla, CA, United States⁵).

The bioinformatical analyses represented by Shel and Python scripts as well as their results are deposited at Zenodo⁶.

Results

Survival Assays Show High Susceptibility to Heat Stress of C. lari and C. coli

The two strains C. coli RM2228 and C. lari RM2100 were heat stressed by cultivation at 46°C under microaerobic conditions and survival of these strains was determined over a period of 24 h. Experiments were started with approximately 7 log₁₀ CFU/ml (Figure 1A). Cell counts of C. coli and C. lari were reduced by one log level after 2 h of incubation at 46°C. C. lari cell counts decreased to 3.4 log₁₀ CFU/ml after 4 h and further 1.1 log₁₀ CFU/ml after 6 h of heat stress. The cell counts of C. coli were reduced to 4.6 log₁₀ CFU/ml after 4 h and further 1.9 log₁₀ CFU/ml after 6 h of cultivation at 46°C. Already after 8 h of heat stress, no surviving cells were observed for both strains. To verify that the observed heat stress survival ability is not only strain-specific, the experiments were repeated with two more field isolates of each species (Figure 1B). Even though variations between the isolates were observed, the overall survival was comparable to the results obtained by the laboratory strains shown in Figure 1A.

FIGURE 1

Figure 1. Survival abilities of C. coli and C. lari at 46°C. Bacterial cultures of approximately 7 log₁₀ CFU/ml were incubated at 46°C and cell counts determined at indicated time points. Data were shown for (A) reference strains and (B) wild type strains of C. coli and C. lari as median ± interquartile range of three biological replicates. Gray: sequenced strains, white: first field strain, and black: second field strain.

Gene Expression Profiling of Known Heat Stress Genes Using qRT-PCR

The expression level of selected genes, known to be involved in the heat shock response of C. jejuni, was investigated over a time period of 60 min of heat stress at 46°C by quantitative RT-PCR in C. coli and C. lari. In both strains increased log₂ fold changes of the expression level for the chaperon genes clpB, grpE, dnaK, groEL, groES, cbpA and the negative transcriptional regulator hrcA were determined, while expression level of dnaJ was only slightly increased in C. coli and not regulated in C. lari (Figure 2). After 15 min of heat stress, the expression level of these genes was highly induced in both strains. While the high expression levels were stable until the end of the experiments (60 min) in C. coli, the levels slightly decreased after 30 min in C. lari. Overall, higher fold changes of expression level of these genes were detected in C. coli compared to C. lari.

FIGURE 2

Figure 2. Gene expression of known heat shock genes in C. coli RM2228 and C. lari RM2100 at 46°C compared to 37°C by qRT-PCR. Gene expression of selected heat shock genes was analyzed by qRT-PCR and normalized to the expression of rpoA and thiC. Fold changes of expression were calculated by the ΔΔCq method of CFX manager software. Log₂ fold changes for each biological replicate were shown and the means connected by the line.

Global Transcriptome Analysis Upon Heat Stress Using RNA-Seq

To get further insight into the heat stress response of both species, whole transcriptome analysis of C. coli and C. lari was performed by RNA sequencing. Therefore, RNA was extracted for each strain after 30 min cultivation at 37 and 46°C of biological duplicates, libraries prepared as described in Section “Materials and Methods” and transcriptome libraries sequenced by Illumina. Log₂ fold changes of expression level were determined by DESeq2 (mapping statistics for all libraries are shown in Supplementary Table S2 and log₂ fold changes of expression level in Supplementary Table S3). The genes already analyzed by qRT-PCR showed similar amplification in the RNA-seq analysis. Expression level of further 18 genes per strain was investigated by qRT-PCR, to correlate the fold changes determined by RNA-seq versus qRT-PCR resulting in a Pearson correlation coefficient R² of 0.6644 for C. coli and 0.5599 for C. lari (Supplementary Figure S1).

After 30 min of heat stress 17.2% of the C. coli (338/1967) and 19.4% of the C. lari genes (300/1545) were shown by RNA-seq analysis to be differentially expressed (Table 1). Of these differentially expressed genes, 67.1% were up-regulated in C. coli, while only 43.6% were up-regulated in C. lari (Figure 3).

TABLE 1

Table 1. Gene expression profile by RNA sequencing of C. coli RM2228 and C. lari RM2100 after 30 min heat stress at 46°C compared to 37°C.

FIGURE 3

Figure 3. Gene expression profile of the whole transcriptome of C. coli RM2228 and C. lari RM2100 after 30 min of heat stress. Whole transcriptome was analysis after 30 min incubation at 46°C by RNA sequencing and gene expression shown in volcano plots as log₂ fold change versus -log₁₀ p-values of biological duplicates. Gray: significantly different expressed genes.

Enrichment of Regulated Genes After Heat Stress in Several Functional Categories

To analyze if genes with similar function were comparably expressed in both species and if a specific class of genes is significantly enriched upon heat stress, a functional categorization based on orthologous groups according to the eggNOG database was used. Nearly half of the differentially expressed genes (55% for C. coli and 49% for C. lari) has not been assigned to any category or their functions are only poorly characterized (Supplementary Table S4). The other differentially expressed genes are distributed over the 20 functional categories, with significant enrichment in the category T (signal transduction mechanisms) for C. coli (27.3%) and significant enrichment in the category I (lipid transport and metabolism) and category Q (secondary metabolites biosynthesis, transport, and catabolism) for C. lari (33.3 and 40%) (Figure 4A).

FIGURE 4

Figure 4. Genes differentially expressed after heat stress according to functional groups. The genes differentially expressed after 30 min of heat stress were shown according to their functional classification (obtained from the evolutionary genealogy of genes: Non-supervised Orthologous Groups http://eggnog.embl.de/version_4.0.beta) for C. coli RM2228 and C. lari RM2100. (A) Total number of regulated and unregulated genes assigned to each functional category were shown for C. lari and C. coli. Significant enrichment for each category were calculated by Fisher’s exact test (*p < 0.05). (B) The figure was created with the open source tool Circos (http://circos.ca). The thickness of each connection displays the log₂ fold change of gene expression whereas the ending of connection displays if the gene expression is up- or down-regulated.

In both species, the majority of differentially expressed genes in category O (posttranslational modification, protein turn-over, and chaperones) and category H (coenzyme transport and metabolism) were up-regulated, while the majority of differentially expressed genes of category J (translation, ribosomal structure and biogenesis) and U (intracellular trafficking, secretion and vesicular transport) were down-regulated (Figure 4B). Expression level of genes belonging to the category of cell wall/membrane/envelop biosynthesis (M), energy production and conversion (C), and amino acid as well as nucleotide transport and metabolism (E, F) were mostly up-regulated in C. coli but down-regulated in C. lari. In contrast, nearly no genes belonging to category cell cycle control, cell division, chromosome partitioning (D) or defense mechanisms (V) were affected in their expression level in both species (Supplementary Table S4).

C. coli and C. lari Share a Low Number of Differentially Expressed Orthologs Genes With Each Other After Heat Stress

To compare the heat stress response of both species in more detail, orthologs between the genomes of C. coli RM2228 and C. lari RM2100 were determined. As the knowledge about gene functionality is more comprehensive for C. jejuni, we also included the genome of C. jejuni NCTC11168 in the ortholog assignment. Altogether 1372 orthologous groups (consisting of two or three genes) could be assigned in the three combined genomes (Supplementary Figure S2) with 759 orthologous genes present in all three genomes. While C. jejuni and C. coli additionally shared 554 orthologous genes, C. lari shared only 18 orthologous genes with C. coli and 41 with C. jejuni, respectively.

For 70.1% of the 338 regulated genes in C. coli, and 57.7% of the 300 regulated genes in C. lari, no ortholog could be determined in the other strain. However, the expression of 35 genes, 29 orthologs, and six genes with the same name, but not defined as orthologs at a BLAST e-value of 1e-6 (marked by asterisks in Table 2), was similarly regulated in both species, while 2 genes showed oppositional expression between the two strains (Table 2).

TABLE 2

Table 2. Similar regulated genes after heat stress in C. coli RM2228 and C. lari RM2100.

Orthologous genes with up-regulated expression levels in both species include the heat stress response related genes hrcA, grpE, dnaK, groEL, groES, clpB, cbpA, hspR, the invasion antigen ciaB, flagellar associated genes flaG and pseA, membrane associated genes (lolA, kefB, Cla_1506/CCO0311), the threonyl-tRNA-synthetase thrS and the exinuclease uvrB.

Orthologous genes down-regulated in both species included translation and ribosomal associated genes (rpmF, rpsL, dusB, aat, queA, an endoribonuclease), the Campylobacter transformation system proteins ctsE and ctsD, the outer membrane efflux protein cmeD, and the flagellar P-ring protein flgI.

The expression of the genes encoding a histidine triad protein (Cla_0423/CCO0504) and a GTP cyclohydrolase (Cla_0006/CCO0022) were down-regulated in C. lari but up-regulated in C. coli (Table 2).

Genes Affected in Expression After Heat Stress in One Species Only

C. coli

Genes assigned to a functional category and with differential expression after heat stress in C. coli only are shown in Table 3 (log₂ fold changes and p-values are shown in Supplementary Table S3). The expression of the transcriptional regulators cmeR, furR1 and a regulator belonging to the ArsR-family were also up-regulated in C. coli, while expression of CCO1284 (a regulator belonging to the Baf-family) and the transcription termination gene nusA were down-regulated. Further, the expression level of several signal transducer or two-component system encoding genes (e.g., cheV, the flgR-ortholog rrp-2, dccS) was up-regulated.

TABLE 3

Table 3. Genes differentially regulated after heat stress in C. coli RM2228 only.

Despite the genes of the translational machinery affected in both species, the expression of the 16S rRNA modification protein (rimM) and the translation termination factor prfA were also down-regulated in C. coli. In contrast, the expression of several ribosomal protein encoding genes as well as the translation initiation factor infC were up-regulated.

Further, increased expression levels for the genes ung and recA, both involved in DNA repair, as well as for the chaperone gene htrA have been determined.

In addition, the expression levels of several genes involved in amino acid (e.g., glnA, bisZ, glnQ, CCO1002), carbohydrate (e.g., fucP, CCO0582 + 83), nucleotide, and coenzyme transport and metabolism were up-regulated after the heat shock in C. coli.

The expression levels of the majority of genes belonging to the cell wall/membrane/envelope biogenesis (e.g., pseF, pseG) and the multidrug efflux system encoding genes cmeA and cmeB were significantly up-regulated whereas the expression of cmeC was only slightly increased (Supplementary Table S3).

The expression level of the flagellar filament structure genes flaA and flaC, and the chemotaxis gene cheW were up-regulated. In contrast, genes belonging to secretory function of the flagella like fliR and fliI as well as the basal body flgG_1 and the energy delivering protein encoded by fliI were down-regulated (Table 3 and Supplementary Table S3).

In contrast, the expression level of other genes involved in energy production and conversion were all up-regulated (Table 3 and Supplementary Table S3).

C. lari

Genes assigned to a functional category and with differential expression after heat stress in C. lari only are shown in Table 4 (log₂ fold changes and p-values are shown in Supplementary Table S3). The expression of the transcriptional regulator genes Cla_1081 and Cla_0872 was up-regulated while rpoN and Cla_0784 (transcriptional regulator of the Crp/Fnr-family) was down-regulated in C. lari (Table 4). Of the two-component system or signal transducer encoding genes, a probable dccS (Cla_1102) and Cla_1204 were up-regulated, while Cla_0781 was down-regulated.

TABLE 4

Table 4. Genes differentially regulated after heat stress in C. lari RM2100 only.

Even though decreased expression levels for the majority of genes belonging to the translational machinery were determined, the expression level of prfB, involved in termination of translation, and Cla_0258, a ribosomal modification protein, were up-regulated (Table 4 and Supplementary Table S3).

The expression of the chaperon genes dnaJ, the DnaJ domain protein encoding gene Cla_0590, as well as xth (involved in DNA repair) was up-regulated. However, several other genes also involved in DNA repair and replication were down-regulated (Table 4 and Supplementary Table S3).

A significant enrichment in the category of lipid transport and metabolism has been determined and the majority of the differentially expressed genes were down-regulated.

The expression of genes involved in LOS biosynthesis (kdsB), rod shape formation (mreC), peptidoglycan biosynthesis (pbpA, murC), or general cell wall/membrane/envelop biogenesis (pglF and Cla_0642) were down-regulated in C. lari.

The flagellin modification genes ptmB and pseB were up-regulated, while the flagellar motor switch protein encoding gene fliN, and the basal body rod protein encoding gene flgG were down-regulated.

Discussion

It has previously been shown that Campylobacter spp. are rather sensitive to elevated temperatures. Baserisalehi et al. (2006) determined similar D-values for C. coli and C. lari in food matrices at elevated temperatures. Our data showed that C. coli and C. lari are similarly susceptible to heat stress at 46°C in Brucella broth, too. In contrast to C. jejuni, studies about the heat shock response of these two Campylobacter species on mRNA or protein level were missing so far. Therefore, we analyzed changes in the whole transcriptome of the two strains after heat stress. The overall percentage of differentially expressed genes was comparable between C. coli (17.2%) and C. lari (19.4%). However, the majority of differentially expressed genes were up-regulated in C. coli (67.1%) but down-regulated in C. lari (56.4%), indicating a different heat stress response of the two species.

Of the overall transcriptional changes detected, only 35 genes were similarly expressed indicating a general role in heat stress response of Campylobacter spp. This group includes the transcriptional regulators hspR and hrcA as well as the chaperons grpE, dnaK, groEL, groES, clpB, and cbpA. These genes are included in the 30 genes whose expression has been most highly up-regulated in both species (Supplementary Table S3), and are also involved in the heat stress response of C. jejuni (Stintzi, 2003; Holmes et al., 2010). Stintzi (2003) described a high and rapid up-regulation of these genes in C. jejuni after temperature increase from 37 to 42°C. This up-regulation of gene expression was only transient, as expression levels nearly reached basal levels after 50 min of temperature up-shift indicating that adaptation to the new conditions was achieved (Stintzi, 2003). As we did not observe any growth at 46°C for the strains investigated in our study, we did not expect adaptation to this conditions, and therefore also no transient changes in the gene expression pattern. This was also supported by the observation that known heat shock genes were still highly expressed after 60 min of heat stress in our study (Figure 2).

Regarding chaperone DnaJ, Konkel et al. (1998) showed higher protein levels after heat stress for C. jejuni. In our study, the expression pattern for dnaJ was contradictory. Slightly increased expression of dnaJ has been determined by RT-PCR for C. coli and by RNA-seq for C. lari only. However, Holmes et al. (2010) suggested that CbpA, showing 44% similarity to DnaJ, is the main DnaJ protein in C. jejuni. In concordance with that, we observed significantly up-regulated expression of cbpA in C. coli and C. lari in both assays. The only gene belonging to the category of posttranslational modification, protein turnover, chaperons down-regulated in both species was aat. The protein Aat might be involved in N-end rule pathway of protein degradation, supposing changes in half-life time of several proteins (Varshavsky, 2011).

The periplasmic chaperone HtrA, which degrades and prevents aggregation of misfolded periplasmic proteins, has been described to be essential for C. jejuni growth at 44°C (Baek et al., 2011). Increased expression of htrA was also observed in C. coli but not in C. lari, suggesting that other chaperons or proteases replace the activity of HtrA in C. lari. Interestingly, also racR was not differentially expressed in C. coli and C. lari, indicating further regulatory differences compared to the heat shock response of C. jejuni (Apel et al., 2012).

Regarding the ribosomal genes, Stintzi (2003) described a transient down-regulated expression in C. jejuni after temperature increase from 37 to 42°C. The expression of these genes decreased until 20 min and nearly reached baseline level at 50 min after temperature increase, suggesting that C. jejuni transiently reshuffled energy for stress damage repair and adaptation to the new growth condition (Stintzi, 2003). Similarly, the majority of differentially expressed genes assigned to the category of translation and ribosomal structure was also down-regulated in C. coli and in C. lari at 30 min after heat stress of 46°C. Also comparable to the observations described for C. jejuni (Stintzi, 2003), the majority of all differentially expressed genes was up-regulated in C. coli, implying high energy costs of heat stress response. According to this hypothesis, the expression of the majority of genes involved in energy metabolism was also up-regulated in C. coli. In contrast, the majority of all differentially expressed genes (including genes assigned to the category of energy metabolism) were down-regulated in C. lari, suggesting a divergent energy saving behavior compared to C. coli under the heat stress condition investigated in our study. The differences in energy metabolism is further supported by the fact that multiple genes involved in energy metabolism and respiration in C. jejuni are missing in C. lari (Miller et al., 2008, 2014). Nevertheless, the question which transcriptional changes would be further transferred to protein level remains unanswered.

Additional regulators seem to be involved in the heat stress response in both strains tested. While C. coli also up-regulated the expression of the transcriptional regulator fur1, the ortholog perR was not regulated in C. lari. PerR has been shown to regulate the expression of at least 104 genes in C. jejuni, most of them involved in oxidative stress response (Palyada et al., 2009). Further, perR expression is up-regulated after acid stress (Reid et al., 2008). Our data could indicate a further role of PerR in the heat stress response of C. coli but not of C. lari.

The two-component system DccRS is required for initial adaptation of C. jejuni to the gastrointestinal milieu of chickens and controls a regulon of six genes annotated as putative periplasmic or membrane proteins, and dccR itself (Wosten et al., 2010). The expression of the putative sensor dccS (CCO1300) was up-regulated after heat stress in C. coli. In concordance, six genes defined as orthologs or annotated as probably C. jejuni genes of the DccRS regulon (CCO0021, CCO0198, CCO0689, CCO0290, CCO1066, CCO1462) were up-regulated in C. coli. For C. jejuni, it has been demonstrated that neither temperature of 42°C nor copper but growth phase is responsible for activation of DccRS (Wosten et al., 2010). However, our data suggest that the DccRS system is also involved in the heat stress response of C. coli at 46°C.

Likewise to C. coli, the expression of further possible transcriptional regulators was regulated in C. lari. The highest up-regulated expression in C. lari was observed for Cla_0805, encoding a peptidase S24 family protein, predicted as transcriptional repressors of the SOS-response in many bacteria (Butala et al., 2009). Interestingly, many of the SOS-response genes described in other bacteria, including the peptidase S24 LexA, are lacking in C. jejuni and no ortholog has been detected in C. coli (Zgur-Bertok, 2013). Whether C. lari has a complete SOS-response has to be elucidated in further studies.

Further changes observed after heat stress include metabolism pathways. Glutamine is the major nitrogen donor for C. jejuni and C. jejuni 81-176, mutated in the glutamine ABC-Transporter (Cj0469), was more susceptible to heat stress compared to the wild type (Lin et al., 2009). Further, C. jejuni use glutamine as a carbon source (van der Hooft et al., 2018). We determined enhanced gene expression of genes involved in glutamine metabolism for C. coli, like the degenerated ammonium transporter (CCO0599), the glutamine synthetase glnA, the glutamine transporter ATPase glnQ, as well as the adjacent located amino acid transporter permease (CCO1002).

Recently, an energy metabolism pathway involved in the formation of pyruvate and L-lactate from L-fucose has been described for C. jejuni (Stahl et al., 2012). These genes are located on a genomic island spanning from Cj0480-89 and determined as ortholog for C. coli RM2228 by our comparison (CCO0578-87) (Muraoka and Zhang, 2011; Stahl et al., 2011). As the expression of these genes was up-regulated, we suggest that the activation of nitrogen assimilation and the L-fucose metabolic pathway is part of heat stress response in C. coli RM2228.

Many genes assigned to the functional category of cell wall/membrane/envelope biogenesis and hypothetical membrane proteins were differentially expressed in C. coli and C. lari. However, the expression of the majority of these genes was up-regulated by C. coli and down-regulated by C. lari suggesting different modifications of cell wall/membrane proteins after heat stress in both species.

Conclusion

Despite similar phenotypic survival abilities at 46°C, distinct variations in the transcriptomic response to heat stress could be determined for C. coli and C. lari. Similar expression changes for chaperone genes after heat stress, which have also been described for C. jejuni, indicate a general response mechanism of Campylobacter species. However, the differences in expression of transcriptional regulators and metabolism related genes observed between the tested C. coli and C. lari strains could be explained in part by the phylogenetic distance and therefore the large differences in the genomic content of these two species. Higher similarities have been determined for the heat stress response of C. coli with already published data for C. jejuni. This is in agreement with the higher amount of orthologous genes determined for C. jejuni and C. coli compared to C. lari.

Overall, the described transcriptomic changes induced after heat stress in C. coli and C. lari indicate huge differences in the heat stress response between the Campylobacter species.

Data Availability Statement

The datasets generated for this study can be found in the GEO Series accession number GSE67486 (www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE67486).

Author Contributions

GG and TA planned the study. CR and CP performed the experiments in this study. KF and CS performed the RNA sequencing and bioinformatic analyses. CR, TA, and GG analyzed results, formatted the data, and drafted the manuscript. TA provided funding. All authors read and approved the final manuscript.

Funding

This work was funded by the Federal Ministry of Education and Research (BMBF) under project number 01Kl1725A (PAC-Campy) as part of the Research Network Zoonotic Infectious Diseases and project number 01Kl1012 (FBI-Zoo).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.00523/full#supplementary-material

FIGURE S1 | Correlation of log₂ fold changes of gene expression by qRT-PCR and RNA-seq. The correlation of the log₂ fold changes of selected gene expression after 30 min of heat stress determined by RNA-seq (x-axis) and qRT-PCR (y-axis) is depicted with corresponding R² values for C. coli RM2228 and C. lari RM2100.

FIGURE S2 | Orthologous mapping of the complete genomes from C. coli RM2228, C. lari RM2100, and C. jejuni NCzTC11168. Orthologous were defined by bidirectional best BLAST-hit search on nucleotide level with max. e-value of 1e-6, word size of 20, and a minimal length of 60%. The Venn diagram shows orthologous gene shared by the C. coli strain RM2228, C. lari strain RM2100, and C. jejuni strain NCTC11168.

TABLE S1 | Primers used in this study. The target gene, primer sequence and concentration, the amplicon length, and annealing temperature of each primer used for qRT-PCR for both strains are indicated in this table.

TABLE S2 | Mapping statistics for all libraries in the RNA-seq. This table indicates the total number of sequenced cDNA reads considered in the analysis and used for alignment, the total number of aligned and of uniquely aligned reads, the total number of alignments, and the percentage of aligned and uniquely aligned reads for the replicates at 37 and 46°C for both strains.

TABLE S3 | Excel file containing the Orthologs of differentially regulated genes after heat stress in C. coli RM2228 and C. lari RM2100. All genes for which a significant regulation of expression was calculated by DESeq2 1.6.3 after 30 min heat stress by RNA-seq analysis were indicated in this table. For each gene, the locus_tag, gene name, gene product, log₂ fold change of gene expression with adjusted p-value, assignment to functional category by eggnog, and their corresponding ortholog are listed. Gene names in parentheses were transferred from orthologs. -: no gene name, function, or category assigned.

TABLE S4 | Percentage of differentially expressed genes after heat stress according to eggNOG functional categories. The table indicates the percentage of regulated genes for each strain and category. Genes assigned in two categories were listed in each category. Significant enrichment was calculated by Fisher’s exact test. ^∗p < 0.05.

TABLE S5 | Strains used in this study. The ID, origin, description, and source of each strain used in this study are indicated in this table.

Footnotes

1. ^ http://frodo.wi.mit.edu/
2. ^ http://hannonlab.cshl.edu/fastx_toolkit/
3. ^ http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE67486
4. ^ http://eggnog.embl.de/version_4.0.beta/downloads.v4.html
5. ^ http://www.graphpad.com
6. ^ https://doi.org/10.5281/zenodo.3669112

References

Alter, T., Bereswill, S., Glunder, G., Haag, L. M., Hanel, I., Heimesaat, M. M., et al. (2011). Campylobacteriosis of man: livestock as reservoir for Campylobacter species. Bundesgesundheitsblatt Gesundheitsforschung Gesundheitsschutz 54, 728–734. doi: 10.1007/s00103-011-1289-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Alter, T., and Scherer, K. (2006). Stress response of Campylobacter spp. and its role in food processing. J. Vet. Med. B Infect. Dis. Vet. Public Health 53, 351–357. doi: 10.1111/j.1439-0450.2006.00983.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Apel, D., Ellermeier, J., Pryjma, M., Dirita, V. J., and Gaynor, E. C. (2012). Characterization of Campylobacter jejuni RacRS reveals roles in the heat shock response, motility, and maintenance of cell length homogeneity. J. Bacteriol. 194, 2342–2354. doi: 10.1128/JB.06041-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Baek, K. T., Vegge, C. S., Skorko-Glonek, J., and Brondsted, L. (2011). Different contributions of HtrA protease and chaperone activities to Campylobacter jejuni stress tolerance and physiology. Appl. Environ. Microbiol. 77, 57–66. doi: 10.1128/AEM.01603-10

PubMed Abstract | CrossRef Full Text | Google Scholar

Baserisalehi, M., Bahador, N., and Kapadnis, B. P. (2006). Effect of heat and food preservatives on survival of thermophilic Campylobacter isolates in food products. Res. J. Microbiol. 1, 512–519. doi: 10.3923/jm.2006.512.519

CrossRef Full Text | Google Scholar

Blomberg, P., Wagner, E. G., and Nordstrom, K. (1990). Control of replication of plasmid R1: the duplex between the antisense RNA, CopA, and its target, CopT, is processed specifically in vivo and in vitro by RNase III. EMBO J. 9, 2331–2340. doi: 10.1002/j.1460-2075.1990.tb07405.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Butala, M., Zgur-Bertok, D., and Busby, S. J. (2009). The bacterial LexA transcriptional repressor. Cell Mol. Life. Sci. 66, 82–93. doi: 10.1007/s00018-008-8378-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Edgar, R., Domrachev, M., and Lash, A. E. (2002). Gene expression omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 30, 207–210. doi: 10.1093/nar/30.1.207

PubMed Abstract | CrossRef Full Text | Google Scholar

EFSA, and ECDC (2014). The European union summary report on trends and sources of zoonoses, zoonotic agents and food-borne outbreaks in 2012. EFSA J. 12:312.

Google Scholar

Forstner, K. U., Vogel, J., and Sharma, C. M. (2014). READemption-a tool for the computational analysis of deep-sequencing-based transcriptome data. Bioinformatics 30, 3421–3423. doi: 10.1093/bioinformatics/btu533

PubMed Abstract | CrossRef Full Text | Google Scholar

Fouts, D. E., Mongodin, E. F., Mandrell, R. E., Miller, W. G., Rasko, D. A., Ravel, J., et al. (2005). Major structural differences and novel potential virulence mechanisms from the genomes of multiple Campylobacter species. PLoS Biol. 3:e15. doi: 10.1371/journal.pbio.0030015

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoffmann, S., Otto, C., Kurtz, S., Sharma, C. M., Khaitovich, P., Vogel, J., et al. (2009). Fast mapping of short sequences with mismatches, insertions and deletions using index structures. PLoS Comput. Biol. 5:e1000502. doi: 10.1371/journal.pcbi.1000502

PubMed Abstract | CrossRef Full Text | Google Scholar

Holmes, C. W., Penn, C. W., and Lund, P. A. (2010). The hrcA and hspR regulons of Campylobacter jejuni. Microbiology 156, 158–166. doi: 10.1099/mic.0.031708-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Hwang, S., Jeon, B., Yun, J., and Ryu, S. (2011). Roles of RpoN in the resistance of Campylobacter jejuni under various stress conditions. BMC Microbiol. 11:207. doi: 10.1186/1471-2180-11-207

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J.-C., Oh, E., Kim, J., and Jeon, B. (2015). Regulation of oxidative stress resistance in Campylobacter jejuni, a microaerophilic foodborne pathogen. Front. Microbiol. 6:751. doi: 10.3389/fmicb.2015.00751

PubMed Abstract | CrossRef Full Text | Google Scholar

Konkel, M. E., Kim, B. J., Klena, J. D., Young, C. R., and Ziprin, R. (1998). Characterization of the thermal stress response of Campylobacter jejuni. Infec Immun. 66, 3666–3672.

PubMed Abstract | Google Scholar

Lin, A. E., Krastel, K., Hobb, R. I., Thompson, S. A., Cvitkovitch, D. G., and Gaynor, E. C. (2009). Atypical roles for Campylobacter jejuni amino acid ATP binding cassette transporter components PaqP and PaqQ in bacterial stress tolerance and pathogen-host cell dynamics. Infect. Immun. 77, 4912–4924. doi: 10.1128/IAI.00571-08

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550.

PubMed Abstract | Google Scholar

Miller, W. G., Wang, G., Binnewies, T. T., and Parker, C. T. (2008). The complete genome sequence and analysis of the human pathogen Campylobacter lari. Foodborne Pathog. Dis. 5, 371–386. doi: 10.1089/fpd.2008.0101

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, W. G., Yee, E., Chapman, M. H., Smith, T. P. L., Bono, J. L., Huynh, S., et al. (2014). comparative genomics of the Campylobacter lari group. Genome Biol. Evol. 6, 3252–3266. doi: 10.1093/gbe/evu249

PubMed Abstract | CrossRef Full Text | Google Scholar

Muraoka, W. T., and Zhang, Q. (2011). Phenotypic and genotypic evidence for L-fucose utilization by Campylobacter jejuni. J. Bacteriol. 193, 1065–1075. doi: 10.1128/JB.01252-10

PubMed Abstract | CrossRef Full Text | Google Scholar

Palyada, K., Sun, Y. Q., Flint, A., Butcher, J., Naikare, H., and Stintzi, A. (2009). Characterization of the oxidative stress stimulon and PerR regulon of Campylobacter jejuni. BMC Genomics 10:481. doi: 10.1186/1471-2164-10-481

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, S. F. (2002). The physiology of Campylobacter species and its relevance to their role as foodborne pathogens. Int. J. Food Microbiol. 74, 177–188. doi: 10.1016/s0168-1605(01)00678-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Parkhill, J., Wren, B. W., Mungall, K., Ketley, J. M., Churcher, C., Basham, D., et al. (2000). The genome sequence of the food-borne pathogen Campylobacter jejuni reveals hypervariable sequences. Nature 403, 665–668. doi: 10.1038/35001088

PubMed Abstract | CrossRef Full Text | Google Scholar

Poropatich, K. O., Walker, C. L., and Black, R. E. (2010). Quantifying the association between Campylobacter infection and Guillain-Barre syndrome: a systematic review. J. Health Popul. Nutr. 28, 545–552.

PubMed Abstract | Google Scholar

Reid, A. N., Pandey, R., Palyada, K., Naikare, H., and Stintzi, A. (2008). Identification of Campylobacter jejuni genes involved in the response to acidic pH and stomach transit. Appl. Environ. Microbiol. 74, 1583–1597. doi: 10.1128/AEM.01507-07

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheppard, S. K., Dallas, J. F., Wilson, D. J., Strachan, N. J. C., McCarthy, N. D., Jolley, K. A., et al. (2010). Evolution of an agriculture-associated disease causing Campylobacter coli Clade: evidence from national surveillance data in scotland. PLoS One 5:e15708. doi: 10.1371/journal.pone.0015708

PubMed Abstract | CrossRef Full Text | Google Scholar

Stahl, M., Butcher, J., and Stintzi, A. (2012). Nutrient acquisition and metabolism by Campylobacter jejuni. Front. Cell Infect. Microbiol. 2:5. doi: 10.3389/fcimb.2012.00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Stahl, M., Friis, L. M., Nothaft, H., Liu, X., Li, J., Szymanski, C. M., et al. (2011). L-fucose utilization provides Campylobacter jejuni with a competitive advantage. Proc. Natl. Acad. Sci. U.S.A. 108, 7194–7199. doi: 10.1073/pnas.1014125108

PubMed Abstract | CrossRef Full Text | Google Scholar

Stintzi, A. (2003). Gene expression profile of Campylobacter jejuni in response to growth temperature variation. J. Bacteriol. 185, 2009–2016. doi: 10.1128/jb.185.6.2009-2016.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Hooft, J. J. J., Alghefari, W., Watson, E., Everest, P., Morton, F. R., Burgess, K. E. V., et al. (2018). Unexpected differential metabolic responses of Campylobacter jejuni to the abundant presence of glutamate and fucose. Metabolomics 14:144. doi: 10.1007/s11306-018-1438-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Varshavsky, A. (2011). The N-end rule pathway and regulation by proteolysis. Protein Sci 20, 1298–1345. doi: 10.1002/pro.666

PubMed Abstract | CrossRef Full Text | Google Scholar

Whiley, H., van den Akker, B., Giglio, S., and Bentham, R. (2013). The role of environmental reservoirs in human campylobacteriosis. Int. J. Environ. Res. Public Health 10, 5886–5907. doi: 10.3390/ijerph10115886

PubMed Abstract | CrossRef Full Text | Google Scholar

Wosten, M. M., van Dijk, L., Parker, C. T., Guilhabert, M. R., van der Meer-Janssen, Y. P., Wagenaar, J. A., et al. (2010). Growth phase-dependent activation of the DccRS regulon of Campylobacter jejuni. J. Bacteriol. 192, 2729–2736. doi: 10.1128/JB.00024-10

PubMed Abstract | CrossRef Full Text | Google Scholar

Yura, T., Kanemori, M., and Morita, T. (2000). “The heat shock response: regulation and function,” in Bacterial Stress Responses, eds G. Storz and R. Hengge-Aronis, (Washington, D.C: Amercian Society for Microbiology), 3–18.

Google Scholar

Zgur-Bertok, D. (2013). DNA damage repair and bacterial pathogens. PLoS Pathog. 9:e1003711. doi: 10.1371/journal.ppat.1003711

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M. J., Xiao, D., Zhao, F., Gu, Y. X., Meng, F. L., He, L. H., et al. (2009). Comparative proteomic analysis of Campylobacter jejuni cultured at 37°C and 42°C. Jpn. J. Infect. Dis. 62, 356–361.

Google Scholar

Zhou, Y., Bu, L., Guo, M., Zhou, C., Wang, Y., Chen, L., et al. (2013). Comprehensive genomic characterization of Campylobacter genus reveals some underlying mechanisms for its genomic diversification. PloS One 8:e70241. doi: 10.1371/journal.pone.0070241

PubMed Abstract | CrossRef Full Text | Google Scholar