The Complete Campylobacter jejuni Transcriptome during Colonization of a Natural Host Determined by RNAseq

Campylobacter jejuni is a major human pathogen and a leading cause of bacterial derived gastroenteritis worldwide. C. jejuni regulates gene expression under various environmental conditions and stresses, indicative of its ability to survive in diverse niches. Despite this ability to highly regulate gene transcription, C. jejuni encodes few transcription factors and its genome lacks many canonical transcriptional regulators. High throughput deep sequencing of mRNA transcripts (termed RNAseq) has been used to study the transcriptome of many different organisms, including C. jejuni; however, this technology has yet to be applied to defining the transcriptome of C. jejuni during in vivo colonization of its natural host, the chicken. In addition to its use in profiling the abundance of annotated genes, RNAseq is a powerful tool for identifying and quantifying, as-of-yet, unknown transcripts including non-coding regulatory RNAs, 5’ untranslated regulatory elements, and anti-sense transcripts. Here we report the complete transcriptome of C. jejuni during colonization of the chicken cecum and in two different in vitro growth phases using strand-specific RNAseq. Through this study, we identified over 250 genes differentially expressed in vivo in addition to numerous putative regulatory RNAs, including trans-acting non-coding RNAs and anti-sense transcripts. These latter potential regulatory elements were not identified in two prior studies using ORF-based microarrays, highlighting the power and value of the RNAseq approach. Our results provide new insights into how C. jejuni responds and adapts to the cecal environment and reveals new functions involved in colonization of its natural host.


Introduction
Campylobacter jejuni is an important human pathogen and a leading cause of bacterial derived gastroenteritis worldwide [1,2]. Human Campylobacter infections are sporadic foodborne diseases commonly associated with foods of animal origin, with C. jejuni being the predominant species associated with human illness [3]. Preventative measures to limit human exposure have focused on identifying genes and loci required for colonization [4][5][6][7][8] as these represent targets for anti-Campylobacter strategies. Although much research has been carried out characterizing regulatory mechanisms in this human pathogen, most of this work has focused on in vitro studies with little knowledge generated towards understanding the mechanisms of global gene expression during colonization.
Regulatory networks in C. jejuni have been characterized principally through in vitro methods including quantitative reverse transcriptase PCR (qRT-PCR), transcriptional reporter strains and microarray studies. These have yielded tremendous amounts of information as to how this organism regulates gene expression; however, many of these studies focused on regulatory mechanisms of specific genes and loci in response to specific substrates and culture conditions in vitro. Only two studies have evaluated the transcriptome of C. jejuni during in vivo colonization of an animal model, including the day-of-hatch chicken model [18] and the rabbit ileal loop model [27]. These studies, utilizing microarray approaches, determined that there are marked differences in gene expression profiles between in vivo and in vitro samples, illustrating the importance for continued research into how C. jejuni regulates gene expression during colonization.
High throughput sequencing of cDNA libraries (RNAseq) has emerged as a powerful approach for mapping transcriptomes and profiling gene expression in diverse bacteria [28][29][30][31][32]. RNAseq has several key advantages over microarray analysis, including 1) the ability to detect and quantify transcripts derived from all regions of the genome, 2) a large dynamic range that affords high sensitivity for low-abundance transcripts, and 3) single nucleotide resolution [33].
Regulatory mechanisms involving small, non-coding RNAs (ncRNAs), or riboswitches have not been well-characterized in C. jejuni. Recently, two in vitro RNAseq studies mapped the transcriptome of C. jejuni and other Campylobacter species [29,30]. One study characterized the regulon of the sigma factor RpoN and correlated transcriptome expression data with protein expression profiles [29]. In the other study, RNAseq was used to map and compare transcriptional start sites, characterize promoter structures and analyze CRISPR RNAs in multiple Campylobacter strains and species [30]. Moreover, both studies identified a wide repertoire of potential non-coding RNA species. These in vitro RNAseq studies have yielded much information into regulatory mechanisms encoded by Campylobacter species; however, this powerful tool has yet to be applied to characterizing the transcriptome of C. jejuni during colonization of its natural host. Although, C. jejuni encodes non-coding RNA species, it does not encode the conserved RNA binding protein Hfq [9,10], a key component of sRNA regulatory mechanisms in other organisms [34][35][36][37]. C. jejuni does, however, encode another conserved RNA binding protein CsrA [9,10] which, in other organisms, plays a role in regulating many processes including motility [38], virulence [39], biofilm formation and central carbon metabolism [40]. The C. jejuni CsrA homologue is involved in the oxidative stress response, biofilm formation, and host cell invasion [41]. Unlike other microbes that use CsrA, C. jejuni does not encode any obvious homologue of the CsrA-antagonizing small RNA csrB; thus the mechanism of CsrA regulation in this pathogen is unknown. Although the regulatory mechanisms of non-coding RNAs have yet to be characterized in C. jejuni, the genome of the closely related organism Helicobacter pylori (which like Campylobacter encodes CsrA but not Hfq), encodes noncoding RNAs as well as riboswitches and antisense transcripts involved in regulatory processes [32,42], suggesting C. jejuni may also regulate gene expression with non-coding RNAs.
Here we report the first complete transcriptome of C. jejuni during in vivo colonization of the chicken cecum using RNAseq. We used a strand-specific, library construction protocol to identify both sense and antisense transcripts. Through this study we identified 272 genes that are differentially expressed in vivo compared to in vitro mid-log and stationary phase grown cultures and identified 51 potential RNA regulators including small, non-coding RNAs (ncRNAs), anti-sense transcripts and probable riboswitches. We also identified the structure of the chick cecal microbiome and identified differences in microbiome structure in C. jejuni colonized chicks compared to mock infected controls. Overall, we identified in vivo expression patterns and possible uncharacterized regulatory mechanisms that provide insight into how C. jejuni is able to adapt to diverse environments.

Results and Discussion
Mapping the in vivo and in vitro transcriptome of C. jejuni using Illumina-based RNAseq C. jejuni gene expression is highly regulated; however, the mechanisms of this regulation are still largely uncharacterized.
To profile the C. jejuni transcriptome in vivo, RNA was isolated from the cecal contents of chicks seven days post colonization with an average of ~7.5 x 10 8 C. jejuni CFU/g cecal contents.
To generate enough starting RNA for library construction (5 µg of DNA-free RNA), RNA isolated from five separate chicks housed in the same brooder, was pooled together. In total, three pools of RNA were isolated from three separate sets of chicks. Campylobacter strand specific, bar-coded cDNA libraries were generated using methods described by Mandlik et al. [31]. Bar-coded cDNA libraries were also generated from three independent cultures of C. jejuni grown in Mueller Hinton Broth (MHB) to mid-log and stationary phase, respectively. The six in vitro libraries were sequenced on a total of two lanes of the Illumina HiSeq2000 platform; each of the three in vivo samples was sequenced in its own lane. The average number of reads obtained for in vitro samples was ~125.4 million reads per sample, with ~87.2% of reads mapping to the C. jejuni 81-176 chromosome and two plasmids, pVir and pTet ( Table  1). Sequencing of the in vivo samples resulted in ~378.8 million reads per sample, in which ~14.6% of reads mapped to the C. jejuni chromosome and plasmids (Table 1). Most reads mapped to open reading frames; however, many reads mapped to intergenic non-coding regions on the chromosome as well as locations anti-sense to open reading frames (ORFs), suggesting that C. jejuni encodes ncRNAs that could be involved in uncharacterized regulatory mechanisms ( Figure  1A).

Structural differences in the chicken microbiome do not correlate to changes in C. jejuni global gene expression during colonization
To determine the relative abundance of C. jejuni in the chick ceca during colonization, we performed 454-pyrosequencing of the 16S rRNA gene isolated from the same C. jejuni colonized chicks used in RNAseq experiments, as well as PBS (mock) infected chicks. DNA based 454-pyrosequencing analysis of the V3-V5 region of the 16S rRNA gene yielded a total of 87,663 sequences, which after trimming and chimera removal by the analysis program mothur [43], yielded 48,571 sequences with minimum read length of 250 base pairs. On average, Campylobacters represented between 5-10% of the total number of sequences identified in the ceca of brooder 1 (10%), 2 (5%) and 3 (5%), (Figure 2A). We identified differences in the structure of the cecal microbiota at the phylum, family and OTU (Operational Taxonomic Unit) level between brooder sets and PBS control chicks (Figures 2 and S1). Previous studies investigating the cecal microbiota of pathogen free and C. jejuni colonized chicks showed a difference in phylotype distribution; however, no changes in the microbiota functional gene content was identified [44]. In this study, chicks colonized with C. jejuni had microbiome structures distinct from the PBS control group ( Figure 2B). The structure of the cecal microbiota of chicks colonized with C. jejuni housed in brooder 1 was statistically different from chicks housed in brooders 2 and 3 ( Figure 2B). Moreover, at the family level, chicks in brooder 1 showed an increase (compared to the other brooders) in the relative abundance of Ruminococcaceae, Leuconostocaceae and Enterococcaceae families and a decrease in members of the Lachnospiraceae family ( Figure 2A). Brooder 1 chicks also had an increase in the relative abundance of Campylobacteraceae (~10%) compared to brooders 2 and 3 (~5%). This poses questions as to whether the indigenous cecal microbiota dictates the level of Campylobacter colonization, or if the level of Campylobacter colonization results in distinct shifts in the microbiota. Although there was a significant difference in the bacterial community structure between brooders, there was no significant difference (avg. R-value of 0.979) in the global transcriptome profiles ( Figure S2, G-I), indicating that changes in the cecal microbiota, observed in this study, did not result in changes in the global gene expression of C. jejuni. This is in agreement with the result that changes in the microbiota of chicks colonized with C. jejuni did not alter functional gene content of the microbiome within the chick cecum [44]. In this study, we observed no changes in transcriptome profiles between animal replicates ( Figure S2 G-I); in contrast, a previous transcriptome study of C. jejuni, using the rabbit illeal loop model (RIL), showed wide variances in transcript abundance between animals (R value of 0.49) [27]. The difference between the two sets of findings could be due to differences in animal models (RIL v. 1-day-old chicks), including host immune response and the host microbiota structure and function, and/or the result of differences in transcriptomic technologies (microarray v. RNAseq).

Differential gene expression during colonization of the chicken cecum and validation of RNAseq
To analyze differential expression of transcripts between in vivo and in vitro samples we utilized the variance analysis package DESeq [45]. A substantial number of genes were significantly differentially expressed (>4-fold, p adj <0.05) in vivo compared to both in vitro conditions. Differential expression profiles of chick samples compared to mid-log samples identified 149 genes that were differentially expressed; 135 transcripts showed increased abundance in vivo while 14 genes had decreased abundance ( Figure 1B, Tables S1 and S2). Comparing DESeq of chick samples to stationary samples, 152 genes were significantly differentially expressed (>4-fold, p adj <0.05); 95 transcripts with increased abundance compared to 57 with decreased abundance ( Figure 1B and Tables S3 and S4). Of the genes with increased abundance in vivo, 29 were increased in comparison to both in vitro conditions ( Figure 1B, Table 2). No transcripts were found to be decreased in abundance in both conditions with a >4-fold cutoff. All genes significantly differentially expressed (>4-fold, p adj <0.05) during in vivo colonization compared to in vitro growth, and those differentially expressed between in vitro growth conditions, are listed in Tables S1-S6. Expression profiles of every gene, under all growth conditions, are listed in Table S9.
To validate expression patterns identified by RNAseq, we performed qRT-PCR on RNA isolated from in vitro and in vivo grown cultures. RNA used in qRT-PCR experiments was isolated independent from samples used in RNAseq library construction. Primers were designed to amplify genes with increased, decreased or unchanged levels of transcript abundance in vivo compared to in vitro samples (Table S8). RNAseq identified katA as one of the highest abundant transcripts during colonization compared to both in vitro growth conditions, and this increase in abundance was confirmed by qRT-PCR ( Figure 3). Furthermore, genes identified as having decreased transcript abundance (Cjj0315 and Cjj0067) or no changes in transcript abundance (cjaC and ccoN) during colonization, compared to in vitro grown cultures, were also confirmed by qRT-PCR, validating the RNAseq expression profiles ( Figure 3). Because differential expression analysis  A. Defining members of the cecal microbiota by relative rank abundance plots (>1%) at the bacterial family level. The average relative abundance for each treatment group (PBS control n=5; brooder 1 n=5; brooder 2 n=5; brooder 3 n=5) is represented in the bar plot. B. Principle coordinates analysis (PCoA) illustrating the community structure relationship between chicks from different treatment groups. This ordination was generated using a Yue and Clayton-based distance matrix representing the relative abundance of OTUs in each community at a 3% OTU definition level. The community of each chick is indicated by a colored symbol (PBS control = blue; brooder 1 = black; brooder 2 = green; brooder 3 = red). All brooder bacterial communities were significantly different from PBS controls (p=0.008 for brooder1 and p=0.007 for brooder 2 and 3; AMOVA). Brooder 1 was significantly different from Brooder 2 and 3 (p=0.003 and p=0.015, respectively; AMOVA).  replicates, of all three growth conditions, showed no statistical differences ( Figure S2).

Insights into the intestinal environment during C. jejuni colonization derived from global transcript regulatory patterns
To identify global regulation patterns occurring in vivo, we graphed DESeq analyses on MA plots (Figure 4). MA plots represent the log 2 of the ratio of abundances of each ORF between the indicated conditions [M], plotted against the average log 2 abundance of that ORF in all conditions [A]. Three functional gene groups were identified as increased in abundance during colonization including iron transport, phosphate transport and oxidative and nitrosative stress responses. Iron has been consistently identified as a critical micronutrient for C. jejuni growth and colonization [46][47][48][49][50]. Increased expression of the heme transport genes chuA and chuB were previously identified in an C. jejuni in vivo microarray transcriptome study [18]. Utilizing RNAseq, we confirmed the increased abundance of chuA and chuB transcripts along with chuC and chuD of the heme transport operon (chuABCD). Additionally, we identified two other iron transporters with increased expression in vivo, including the outer membrane TonB-dependent energy transduction system exbB-exbD-tonB and the ferric transport system FTR1-p19. The increased abundance of three iron transport systems suggests that the chicken cecum is an extremely iron-limiting environment.
Transcripts from the phosphate transport system (pstSCAB) were increased in abundance in vivo compared to both in vitro conditions ( Figure 5B, Tables S1 and S3); however, in contrast to the iron transport systems noted above, increased expression of this locus was not observed in either previous in vivo microarray transcriptome study [18,27]. In addition to this ABC-type transport system, C. jejuni encodes a low affinity non-specific phosphate transport system pitAB (Cjj1208, Cjj1209) [51], which was not differentially expressed to any significant extent in our studies. In response to phosphate limitation, the pst transport system is activated by the twocomponent system PhoSR; pst is required for growth in phosphate-limited media [52] while the pitAB locus has yet to be characterized. A role for pitAB and pstSCAB in vivo has yet to be determined in C. jejuni, but increased abundance of the pst system in vivo suggests that the chick cecum is limited in phosphate and that this system could be required for colonization based on its importance in vitro [52].
Furthermore, increased in vivo transcript abundance (compared to in vitro growth) was observed from genes involved in oxidative and nitrosative stress responses ( Figure  4, Figure 5, Tables S1 and S3). The second highest differentially expressed gene in vivo, katA (encoding catalase), is critical for the C. jejuni response to H 2 O 2 and is required for colonization of the chicken cecum ( Figure 5A) [53][54][55]. Additionally, the truncated flavohemoprotein (cgb), involved in the nitrosative stress response [56], was highly expressed in vivo ( Figure 5C). The high level of expression of these two important stress response genes during in vivo colonization corresponds with studies demonstrating that C. jejuni induces a  host immune response in the chicken [57][58][59]. Specifically, RNAseq analysis of chicken cecal tissue during colonization with C. jejuni showed increased expression of NOXO1, (NADPH oxidase organizer 1) which regulates respiratory burst [60]. Moreover, infection of human HCT-8 intestinal epithelial cells with C. jejuni activated host Nox1 and Duox2, resulting in increased H 2 O 2 production and consequent inactivation of CjtK. CjtK, a tyrosine kinase encoded by C. jejuni, activates GalE, required for glycosylation and capsule formation, important factors in adherence and invasion [61]. The high level of catalase gene expression during in vivo growth suggests that C. jejuni is responding to host H 2 O 2 production, which might be a mechanism to avoid inactivation of CjtK allowing for adherence and invasion. Furthermore, C. jejuni elicits a TLR response in chicks, specifically activating TLR-2, TLR-4 and TLR-21 [62]. Most importantly, TLR signaling was shown to be induced when using C. jejuni cell lysates compared to whole cells, proposing that bacterial lysis (evoking a hostile host environment) is required for full TLR activation [62]. Although C. jejuni is considered a commensal that does not induce human-like disease symptoms in chicks, it does elicit an immune response; therefore, increasing expression of stress response genes could possibly result in a mechanism for C. jejuni to postpone lysis, preventing a robust immune response, resulting in persistent colonization. Overall, we identified 272 transcripts that were differentially expressed in vivo (Tables S1-S4). Many of the genes whose abundance increased in vivo contribute to cellular processes including central carbon metabolism, electron transport, biosynthetic processes and transport systems, indicative of a nutrient-restricted environment in the chick cecum compared to the rich in vitro growth medium MHB. We observed that in vivo abundance of RNA transcripts from the pTet plasmid were extremely low, with some transcripts never detected (Figure 4). This suggests either that C. jejuni strongly down-regulates expression of this plasmid in vivo, or that the plasmid is lost at a high rate during colonization. Both pTet and pVir are conjugative plasmids in C. jejuni 81-176 and both contain type-4 secretion system (T4SS) genes, with only genes encoded on pVir functionally characterized [63][64][65]. Involvement of pVir in attachment, invasion and disease during colonization of ferrets has been demonstrated [65]; however, no mechanism by which pTet might contribute to colonization has been uncovered. Genes involved in histidine biosynthesis and nitrogen metabolism (hisH and hisF) [66] were decreased in vivo (Table S2), suggesting that the chick cecum could be replete with histidine and/or glutamate. In C. jejuni, hisH and hisF have been implicated in the biosynthesis of intermediates of the flagellar glyscosylation pathway [67]; however, no role in chick colonization has been determined. Similarly, the major antigenic peptide (PEB3) was also determined to be decreased in abundance in vivo (Table S2). This glycoprotein is immunoreactive [68,69] suggesting that decreased expression could be a mechanism to evade host recognition leading to persistent colonization.
Comparison of RNAseq expression profiles to those of a previous in vivo transcriptome study using microarrays revealed many differences in expression patterns along with some similarities [18]. Most evident are the total number of differentially regulated genes between the two experiments, with microarray identifying only 68 differentially expressed genes (>2-fold), while RNAseq identified 272 differentially expressed genes, even using a higher stringency cutoff (>4fold). Functional groups identified as differentially expressed in both experiments included components of electron transport, central carbon metabolism and transport systems. C. jejuni encodes a highly branched electron transport chain allowing for adaptation to specific environments [7,8,70]. C. jejuni is considered a microaerophile but has been shown to grow anaerobically in the presence of alternative electron acceptors including DMSO and nitrate [8]. RNAseq identified increased abundance of many transcripts associated with electron transport, including DMSO reductase, sulphite oxidase, the low-affinity oxidase (cioAB) and cytochrome c biogenesis (Tables S1-S4). Alternatively, there was decreased in vivo abundance of NADH: ubiquinone oxidoreductase transcripts (~2-3-fold), which encode proteins responsible for using reduced flavodixin as an electron donor instead of NADH [71], and there was no change in the abundance of the high-affinity oxidase CcoN ( Figure 3B). Taken together, and as previously hypothesized, C. jejuni most likely encounters a limited oxygen environment in the cecum, as CioAB is a low-affinity oxidase [72], and C. jejuni has been shown to be able to respire off substrates such as sulphite and DMSO [8,70,73]. We did not observe an increase in expression of the nitrate reductase genes (napAB) as identified in the previous chick in vivo microarray study [18]. This difference in expression of napAB could be a result of growth temperature variation of in vitro grown cultures. The nitrate reductase locus is up-regulated at 42°C, compared to 37°C [74]. In vitro cultures for RNAseq experiments, reported here, were grown at 42°C in an effort to identify genes that were differentially regulated based on host colonization factors and the host nutritional environment, not differences in temperature, as the internal temperature of the chicken is 42°C. Additionally, expression differences could be due to different colonization conditions (including strain differences, inoculum load, duration of colonization and in vitro growth conditions). Finally, these differences could also be attributed to differences in transcriptome mapping technologies, as RNAseq affords direct quantification of transcripts and high sensitivity for low-abundance transcripts. As transcriptome studies are snapshots of gene expression, it would be interesting to determine C. jejuni gene expression profiles over time during the course of colonization. This would allow for a more in-depth profile of the transcriptome during initial colonization through persistent colonization.

Identification of non-coding and anti-sense RNA species
We next mined our RNAseq data for non-coding RNAs and identified 51 putative ncRNA candidates on the chromosome (Table 3) and on plasmid pTet (Table S7), including transcripts derived from intergenic regions and anti-sense to open reading frames ( Figure 6). Our analysis identified nearly all the C. jejuni ncRNAs annotated in the Rfam database [75], including the thiamin pyrophosphate (TPP) riboswitch, the bacterial signal recognition particle (SRP), tmRNA and the RNA component of RNase P [9,76]. Similarly, we detected expression of 17 of the 25 ncRNAs identified in previous C. jejuni RNAseq studies [29,30]. In several cases, expression of ncRNAs previously identified in C. jejuni strain NCTC11168 were not detected in our analysis of strain 81-176, despite the conservation of these loci in both strains, an observation consistent with recent studies in Campylobacter sp. [30] and Listeria sp. [77] showing that expression of conserved ncRNAs may diverge significantly even among strains of the same species. Additionally, several ncRNAs previously identified in C. jejuni 81-176 [30] were not identified in our study, likely reflecting differences in growth temperatures, RNA library construction protocols, and/or ncRNA prediction algorithms.
The abundance of most of these candidate ncRNAs increased during stationary phase compared to mid-log phase ( Figure 7A), suggesting they are likely induced in response to increased stress and/or nutrient limitation. The abundance of  these putative ncRNAs was also induced in vivo compared to mid-log cultures, though not as much as in stationary phase ( Figure 7B and 7C), again suggesting their expression is likely regulated by similar stimuli and conditions encountered by C. jejuni during its transition from log to stationary growth in culture.
To investigate possible regulatory mechanisms of the newly identified ncRNAs we utilized the sRNA target prediction program TargetRNA [78]. Non-coding RNAs are a unique class of regulators that can regulate gene expression by different mechanisms, including altering mRNA interactions with the ribosome [79] and modulating mRNA stability [80][81][82][83]. Additionally, ncRNAs can act as cis-regulatory elements or as trans-regulatory elements [84]. We chose to focus on the intergenic ncRNAs for TargetRNA analysis, as these are more likely to act as trans-regulatory elements and potentially have multiple targets. TargetRNA identified at least 4 targets (pvalue < 0.01) for eight intergenic encoded ncRNAs ( Table 4). The pattern of abundance of many of these targets was the same or opposite of their cognate ncRNAs among the conditions tested (Table 4), suggesting their stability may be regulated by the ncRNA and lending further credence to the TargetRNA predictions. Those targets whose expression patterns were not well correlated with their cognate ncRNA may represent false predictions or correspond to ncRNA: mRNA interactions that do not significantly alter message stability. Further characterization of these putative ncRNAs and target mRNAs will likely yield key insights into C. jejuni regulatory mechanisms and pathways.
The small size of its genome, its lack of an Hfq homologue, and the relatively limited success of genomics-based approaches for identifying ncRNAs in C. jejuni [85] have led to speculations that C. jejuni does not encode the sizeable arsenal of ncRNAs found in many other bacteria. However, the discovery of numerous putative C. jejuni ncRNAs in this study, and in other recent RNAseq-based studies, suggests that C. jejuni does indeed possess a robust RNA-mediated regulation system and provides insights into the mechanisms by which this organism, which encodes so few canonical transcription factors, is able to effectively adapt to rapidly changing environmental conditions and niches.

The Campylobacter jejuni RNAseq transcriptome browser
The RNAseq data generated during the course of this study contains a wealth of information that will be a valuable resource to the numerous researchers studying C. jejuni. To make these data readily accessible to this community, we have generated coverage plots for each RNAseq dataset that can be visualized in the GenomeView browser [86] using the following link: http://www.broadinstitute.org/software/genomeview/ supplemental/CjVJD13/. This browser allows visualization of expression profiles of all annotated ORFs (under all growth conditions), as well as non-coding RNA species; including 5'-UTRs, anti-sense transcripts and intergenic ncRNAs.

Summary
This is the first report of the complete transcriptome of C. jejuni during colonization of the chicken cecum using RNAseq. Through our strand specific library preparation, we identified 51 probable non-coding RNA species, 29 of which were not previously annotated. We also determined the expression profile of C. jejuni during colonization, identifying 272 genes that are significantly differentially expressed in vivo. Additionally, we characterized the microbiota of the chick cecum during C. jejuni colonization. It was shown that C. jejuni represents ~5-10% of the bacterial population and that structural differences in microbiota profiles do not result in changes in the global C. jejuni transcriptome. C. jejuni must rapidly regulate gene expression in response to diverse environmental conditions; however, the regulatory mechanisms underpinning these responses are poorly understood. This work provides key insights into the genes and functions involved in the transition of C. jejuni from logarithmic to stationary growth and during its adaption to the chicken cecal environment, as well as revealing a number of putative regulatory ncRNAs that may play key roles in these adaptations.

Bacterial strains and growth conditions
C. jejuni strain DRH212 (an 81-176 streptomycin resistant derivative [87]) was grown on Mueller Hinton agar (BD, Sparks, MD) supplemented with 10% sheep blood (BA) or in Mueller Hinton Broth (MHB) (BD, Sparks, MD). C. jejuni was cultured at 42°C microaerobically in a tri-gas incubator, constantly maintained at 5% O 2 , 10% CO 2 , balanced with N 2 . Liquid cultures of C. jejuni were grown in MHB to mid-log or stationary phase (as monitored by optical density Ab 600 ).

Chicken Colonization
Animal work carried out in this study followed the recommendations of the National Institutes of Health in the Guide for the Care and Use of Laboratory Animals. This protocol was approved by the University Committee on Use and Care of Animals at the University of Michigan (Protocol 10462-1). All efforts were made to minimize suffering of animals throughout the course of the experiment.
One-day-old chicks were inoculated orally with 1x10 3 CFU of C. jejuni in PBS. Cells used in the inoculum were grown on BA for ~16 h, washed and resuspended in PBS to 10 4 CFU/mL. Chicks were inoculated with 100 µl of cell suspension via oral gavage. Seven days post inoculation, chicks were sacrificed and cecal contents were removed. For RNA isolation, cecal contents were immediately suspended in 1 mL of RNAlater (Qiagen, Valencia CA) and frozen at -80°C. For colonization load and DNA extraction for 454-pyrosequencing, a subsample of cecal contents was removed and diluted in PBS.

RNA isolation
In vitro RNA was isolated from C. jejuni grown in MHB under microaerobic conditions to mid-log and stationary phase (as monitored by optical density Ab 600 ). An equal volume of RNA stop solution (95% EtOH, 5% phenol) was added to the cultures immediately after removal from microaerobic conditions and cells were harvested by centrifugation at 4°C. RNA was extracted with TRIzol® (Invitrogen, Grand Island, NY) The In Vivo Transcriptome of C. jejuni by RNAseq as per manufacturer's instructions. Contaminating DNA was removed with two treatments of TURBO™ DNase (Invitrogen) per manufacturer's instructions. DNase treated RNA was cleaned with RNA clean and concentrate 25 columns after each DNase treatment step (Zymo Research, Irvine, CA). Conformation of DNA removal was assessed by PCR with three sets of primers ranging in amplification products of 100, 150, 250 bp and the Qubit dsDNA high sensitivity assay kit (Life Technologies) (data not shown).
For in vivo RNA isolation, cecal contents were harvested from individual birds and immediately submerged in 1 ml of RNAlater/cecum. C. jejuni was enriched from cecal contents as performed by Jerome et al. in an effort to deplete contaminating eukaryotic cells and cecal debris [88]. After C. jejuni enrichment, RNA was extracted with TRIzol® and contaminating DNA was removed with TURBO™ DNase as described previously. Purified DNA-free RNA from 5 birds (housed in the same brooder) was pooled to create enough starting DNA-free RNA for Illumina cDNA library construction.

Illumina Library Construction
Illumina cDNA libraries were constructed in a similar manner to Mandlik et al. [31]. Five micrograms of RNA was depleted of 16S and 23S rRNA species using the Gram Negative Ribo-Zero™ rRNA Removal Kit (Epicentre, Madison, WI). Removal of contaminating RNA species was assessed with the Agilent Bioanalyzer RNA 6000 nano chip (Agilent Technologies, Santa Clara, CA). Depleted mRNA was fragmented into 100-500 bp species with fragmentation buffer from the GeneChip® clean up module kit (Affymetrix, Santa Clara, CA). First-strand cDNA was synthesized with random hexamers, Actinomycin D and SuperScript III (Life Technologies, Grand Island, NY). Second strand cDNA was synthesized with dUTP replacing dTTP as described by Levin et al [89]. Double stranded cDNA ends were repaired and adenylated as described in the Illumina Truseq™ RNA sample preparation low throughput (LT) protocol (Illumina, San Diego, CA). Bar-coded Illumina adapters were ligated to the ends of the cDNA libraries and adapter-cDNA libraries were treated with Uracil-N-glycosylase (UNG) for 15 min at 37°C, followed by 95°C for 5 min. UNGtreaded cDNA was enriched by low-cycle PCR (8-cylces) with Illumina adapter specific primers. Final cDNA libraries were cleaned with two treatments of AMPureXP beads (Beckman Coulter, Brea, CA), and sequenced using 50 bp paired-end reads on an Illumina HiSeq2000 platform at the University of Michigan Sequencing Core.

Quantitative reverse transcriptase polymerase chain reaction
PCR primers used in this study are listed in Table S8 and were used to amplify nucleotide fragments of genes of interest between 100 and 150 bp. Isolation of RNA from in vitro and in vivo cultures (independent from RNAseq library preparations) was performed as mentioned previously. qRT-PCR was performed by using the Quantitect Sybr, Green RT-PCR kit (Qiagen, Valencia, CA) as by Taveirne et al. [17] and analyzed using the 2(-Delta Delta C(T)) method [90].

RNAseq analysis
Reads were aligned to the C. jejuni chromosome, pTet, and pVir sequence files (RefSeq NC_008787, NC_008770, and NC_008790) using BWA [91] version 5.9. Gene annotations were obtained from RefSeq and Rfam [75]. The overall fragment coverage of genomic regions corresponding to features such as ORFs and rRNAs was conducted as described [31]. Differential expression analysis was conducted using DESeq [45].  [43]. Mothur was used to group or assign bacterial 16S rRNA gene sequences into Operational Taxonomic Units (OTUs) using a 3% species level definition. Classifications were determined by comparing sequences to the Ribosomal Database Project [92]. Classified OTUs were used to determine the relative abundance of bacterial phyla and family in each sample. Principal coordinates analysis (PCoA) was used to assess community similarity among all samples and was calculated based on Yue and Clayton-based distance matrix [93]. An AMOVA (analysis of molecular variance) was performed to determine statistical differences between the community structures of each treatment group. Heatmaps of microbiome data (OTUs) were made using the heatmap command in R [94]. The heatmap command re-orders the data in the rows and the columns separately, so that similar data are grouped together by hierarchical clustering, as shown by dendrograms for the rows and the columns. The Euclidean distance was used in the clustering [95]. Figure S1.

Supporting Information
Structural changes to the chicken cecal microbiota. A heatmap of the top 25 (>1%) operational taxonomic units (OTUs) found in the chick cecal microbiome shows differences between the four-treatment groups (PBS control n=5; brooder 1 n=5; brooder 2 n=5; brooder 3 n=5). The distance between treatment groups was measured by a dendrogram in R statistical program. The PBS controls (blue bars) group very distinctly from chicks colonized with C. jejuni in brooder 1 (black bars) and brooders 2 and 3 (green and red bars). Bacterial phylum, family and OTU number are listed on the heatmap. The heatmap scale ranges from 0 to 100% relative OTU abundance. (TIF)  Table S1.
Genes increased in abundance in vivo compared to in vitro mid-log phase cultures. Listed are genes with increased abundance during in vivo colonization compared to in vitro mid-exponential phase broth grown cultures, as determined by DESeq analysis (materials and methods). Only genes significantly differentially regulated (>4fold difference in abundance, p adj <0.05) are listed. p adj <0.05, is a corrected p-value analogous to a false detection rate of < 5%. Genes are grouped by functional classification and by their C. jejuni 81-176 locus numbers and gene name or function. (DOCX) Table S2.
Genes decreased in abundance in vivo compared to in vitro mid-log phase cultures. Listed are genes with decreased abundance during in vivo colonization compared to in vitro mid-exponential phase broth grown cultures, as determined by DESeq analysis (materials and methods). Only genes significantly differentially regulated (>4fold difference in abundance, p adj <0.05) are listed. p adj <0.05, is a corrected p-value analogous to a false detection rate of < 5%. Genes are grouped by functional classification and by their C. jejuni 81-176 locus numbers and gene name or function. (DOCX) Table S3.
Genes increased in abundance in vivo compared to in vitro stationary phase cultures. Listed are genes with increased abundance during in vivo colonization compared to in vitro stationary phase broth grown cultures, as determined by DESeq analysis (materials and methods). Only genes significantly differentially regulated (>4-fold difference in abundance, p adj <0.05) are listed. p adj <0.05, is a corrected pvalue analogous to a false detection rate of < 5%. Genes are grouped by functional classification and by their C. jejuni 81-176 locus numbers and gene name or function. (DOCX )   Table S4.
Genes decreased in abundance in vivo compared to in vitro stationary phase cultures. Listed are genes with decreased abundance during in vivo colonization compared to in vitro stationary phase broth grown cultures, as determined by DESeq analysis (materials and methods). Only genes significantly differentially regulated (>4-fold difference in abundance, p adj <0.05) are listed. p adj <0.05, is a corrected pvalue analogous to a false detection rate of < 5%. Genes are grouped by functional classification and by their C. jejuni 81-176 locus numbers and gene name or function. (DOCX) Table S5. Genes increased in abundance in vitro mid-log compared to in vitro stationary phase cultures. Listed are genes with increased abundance during in vitro midexponential phase broth grown cultures compared to in vitro stationary phase broth grown cultures, as determined by DESeq analysis (materials and methods). Only genes significantly differentially regulated (>4-fold difference in abundance, p adj <0.05) are listed. p adj <0.05, is a corrected pvalue analogous to a false detection rate of < 5%. Genes are grouped by functional classification and by their C. jejuni 81-176 locus numbers and gene name or function. (DOCX) Table S6. Genes decreased in abundance in vitro mid-log compared to in vitro stationary phase cultures. Listed are genes with decreased abundance during in vitro midexponential phase broth grown cultures compared to in vitro stationary phase broth grown cultures, as determined by DESeq analysis (materials and methods). Only genes significantly differentially regulated (>4-fold difference in abundance, p adj <0.05) are listed. p adj <0.05, is a corrected pvalue analogous to a false detection rate of < 5%. Genes are grouped by functional classification and by their C. jejuni 81-176 locus numbers and gene name or function. (DOCX)  )   Table S9. Results of DEseq analysis of RNAseq data for all annotated C. jejuni ORFs. Each annotated gene in the C. jejuni genome, and two plasmids, is listed in locus number order. DESeq analysis for each gene under each growth condition comparison is given, including the base mean reads for each growth condition and the fold change (and log 2 fold change) between each condition. p adj is the adjusted P values for each gene and resVarA/B are the variance values within each set of biological replicates. (XLS)