Pseudomonas putida KT2440 genome ... - Wiley Online Library

Environmental Microbiology (2011) 13(5), 1309–1326

doi:10.1111/j.1462-2920.2011.02430.x

Pseudomonas putida KT2440 genome update by cDNA sequencing and microarray transcriptomics emi_2430

Sarah Frank,1* Jens Klockgether,1 Petra Hagendorf,2 Robert Geffers,2 Ulrike Schöck,3 Thomas Pohl,3 Colin F. Davenport1 and Burkhard Tümmler1 1 Klinische Forschergruppe, OE 6710, Medizinische Hochschule Hannover, Carl-Neuberg-Str. 1, D-30625 Hannover, Germany. 2 Helmholtz Zentrum für Infektionsforschung, Inhoffenstr. 7, D-38124 Braunschweig, Germany. 3 GATC Biotech AG, Jakob-Stadler-Platz 7, D-78467 Konstanz, Germany. Summary Pseudomonas putida KT2440 is a completely sequenced biosafety strain that has retained its capability to survive and function in the environment. The global mRNA expression profiles of the KT2440 strain grown at 10°C and 30°C were determined by deep cDNA sequencing to refine the genome annotation. Transcriptome sequencing identified 36 yet unknown small non-coding RNAs, 143 novel ORFs in 106 intergenic regions, 42 unclassified genes and eight highly expressed leaderless mRNA transcripts. The genome coordinates of eight genes and the organization of 57 operons were corrected. No overrepresented sequence motifs were detected in the 5⬘-untranslated regions. The 50 most highly expressed genes made up 60% of the total mRNA pool. Comparison of cDNA sequencing, Affymetrix and Progenika microarray data from the same mRNA preparation revealed a higher sensitivity and specificity of cDNA sequencing, a relatively poor correlation between the normalized cDNA reads and microarray signal intensities, and a systematic signal-dependent bias of microarrays in the detection of differentially regulated genes. The study demonstrates the power of next-generation cDNA sequencing for the quantitation of mRNA transcripts and the refinement of bacterial genome annotation. Introduction Pseudomonas putida strain KT2440 (Bagdasarian et al., 1981; Regenhardt et al., 2002) is the best characterized Received 26 August, 2010; accepted 6 January, 2011. *For correspondence. E-mail [email protected]; Tel. (+49) 511 5326721; Fax (+49) 511 5326723.

© 2011 Society for Applied Microbiology and Blackwell Publishing Ltd

1309..1326

saprophytic Pseudomonad that has retained its ability to survive and function in the environment. Strain KT2440 is the first Gram-negative soil bacterium to be certified as a safety strain by the Recombinant DNA Advisory Committee (Federal Register, 1982) and is the preferred host for cloning and gene expression for Gram-negative soil bacteria (Ramos et al., 1987). The KT2440 genome sequence was published in this Journal in 2002 (Nelson et al., 2002). Based on the genome annotation two microarrays were designed by Affymetrix and Progenika for genome-wide gene expression studies (http:// www.affymetrix.com; Yuste et al., 2006). Next-generation sequencing (Metzker, 2010) has recently emerged as an attractive alternative for transcriptome studies (Morozova et al., 2009; van Vliet, 2010). Deep cDNA sequencing by pyrosequencing or sequencing-by-synthesis has been applied to resolve the transcriptome of several microbial species, namely Burkholderia cenocepacia (Yoder-Himes et al., 2009), Helicobacter pylori (Sharma et al., 2010), Listeria monocytogenes (Oliver et al., 2009), Pseudomonas syringae (Filiatrault et al., 2010) and Sulfolobus solfataricus (Wurtzel et al., 2009). In this study we exploited next-generation Illumina sequencing technology (Wang et al., 2009) to refine the current annotation of the KT2440 genome. Transcriptome sequencing data of KT2440 grown at 10°C and 30°C were queried for yet undescribed small RNAs and ORFs and employed to validate predicted operons and gene coordinates. Moreover, we compared the sensitivity and specificity of transcriptome data by taking the same RNA preparations for cDNA sequencing and hybridization of Progenika and Affymetrix microarrays.

Results and discussion P. putida KT2440 growth experiments The KT2440 strain was grown in 1.5 l M9 minimal medium with succinate as carbon source at 30°C in a fermenter until mid-exponential phase. After collection of three 10 ml samples, the culture was cooled down to 10°C within 45 min and 2 h later further samples were harvested. Samples were then processed for RNA preparation, cDNA synthesis and labelling for microarray hybridization as described in the Experimental procedures section.

1310 S. Frank et al. KT2440 deep cDNA sequencing The primary Genome Analyzer II cDNA sequencing data (15–20 million reads per run) were aligned by the ELAND and Bowtie software packages to the KT2440 genome and visualized with the Integrated Genome Browser (Nicol et al., 2009, http://igb.bioviz.org). Figure 1 shows examples of a highly expressed small RNA (sRNA) (crcZ, panel A) (Sonnleitner et al., 2009), a redefined operon (PP0768–PP0772, panel B) and a novel ORF (PP0363.1, panel C). Table 1 summarizes the update of the current annotation of the KT2440 genome in terms of the number

Table 1. Summary of P. putida KT2440 genome update.

Gene category

Number of genes/ operons

New sRNAs New protein coding ORFs New unclassified genes Corrected operons Verified operons

36 143 42 57 154

of yet undescribed sRNAs, ORFs, unclassified transcripts and the number of refined genes and operons.

Small RNAs

Fig. 1. Integrated Genome Browser (IGB) interface showing Illumina cDNA reads (Ganesan et al., 2008). (A) sRNA CrcZ at 30°C, (B) verified operon organisation PP0768-PP0772 at 10°C and (C) the newly identified protein coding ORF PP0363.1 at 10°C.

Based on sequence comparisons with validated sRNAs in other Pseudomonas species (Livny et al., 2006; González et al., 2008; Sonnleitner et al., 2008; 2009; Filiatrault et al., 2010), 22 sRNAs could be unequivocally identified in the KT2440 genome (Table 2a). Besides the tmRNA quality assurance system (Dulebohn et al., 2007), the three sRNAs P26, CrcZ and Psr2 were abundant in P. putida, the latter two of which were downregulated at 10°C. The function of P26 and Psr2 are as yet unknown, but CrcZ has recently been demonstrated in Pseudomonas aeruginosa to be vital for carbon repression control (Sonnleitner et al., 2009). Another 56 abundant transcripts were detected in the intergenic regions, none of which showed homology to any protein in the available databases suggesting that the majority of these transcripts represent small RNAs. Fourteen of the 56 transcripts contain one to five sequence motifs of at least 18 base pairs in length (median: 28 bp; inner quartiles: 26–47 bp; range: 18–756 bp) that occur 2 to 19 times (median: 4) elsewhere in the genome (Table 2b). These sequence repeats could mediate a thermodynamically stable binding of the sRNAs to their potential mRNA targets. The sRNA in IGR2438, for example, contains direct sequence repeats also present at the 3′ end of genes that encode regulatory proteins including the sensor kinase PhoQ (Gooderham and Hancock, 2009; Yan et al., 2009). Thus a common action of the sRNA on the potential targets seems reasonable. Non-random distributions of the potential target sequences were seen for the sRNAs in IGR4451 and IGR4535. Direct repeats of IGR4451 map at the 5′ and 3′ ends of PP3962 and at a third site in IGR3684, while for IGR4535 two sequence motifs are duplicated in PP4739 and PP2508. To resolve the functional implications of the binding of the sRNAs to their targets at the DNA and/or RNA level, thorough wetlab experimentation will be required as demonstrated for example in Pseudomonas fluorescens for RsmY and RsmZ (Lapouge et al., 2007).

© 2011 Society for Applied Microbiology and Blackwell Publishing Ltd, Environmental Microbiology, 13, 1309–1326

Spot42-like (spf) RsmY P26 FMN riboswitch YybP-YkoY T44a RsmZ RgsA/P16 P15 Psr2 P6 PrrF CrcZ S15a P31 P32 SsrA tmRNA TPP riboswitch P1 6S RNA RnpB/P28b Cobalaminb

4

IGR 0124 IGR 0371 IGR 0447 IGR 0531 IGR 0761 IGR 1591 IGR 1625 IGR 1968 IGR 3081 IGR 3541 IGR 4487 IGR 4686 IGR 4697 IGR 4709 IGR 4725 IGR 4725 IGR 4739 IGR 4923 IGR 5184 IGR 5203 PP1327 PP1671

Location 106.24 24.21 3891.26 71.95 66.60 340.76 64.33 119.44 58.11 1889.08 53.48 98.98 3110.03 1199.16 186.28 186.28 7274.98 67.93 37.39 99.84 580.14 21.46

30°C 158.65 19.89 3728.81 47.28 43.76 178.11 137.38 158.47 21.21 121.48 76.87 97.56 399.96 668.35 72.00 72.00 7220.54 6.97 174.03 236.91 946.00 70.49

10°C 1.49 1.17 -1.04 -1.52 -1.52 -1.91 2.14 1.33 -2.74 -15.56 1.44 -1.01 -7.78 -1.79 -2.59 -2.59 -1.01 -9.74 4.65 2.37 1.63 3.28

Fold-change

5 390 008 5 596 334 5 912 448 5 934 663 1 513 067 1 866 934

5 354 623

5 325 379

2 229 914 3 466 159

450 801 537 405 616 513 876 117 1 785 128

Start

5 390 370 5 596 228 5 912 629 5 934 842 1 512 692 1 867 155

5 354 527

5 325 527

2 229 710 3 466 041

450 927 537 470 616 380 875 976 1 785 222

Stop

RFAM predictionc

363 107 182 180 376 222

97

149

195 119

127 66 134 142 95

Length + + + + + + + + + + + + + +

Strand

130 354 450 782 537 397 616 551 876 059 1 785 103 1 822 044 2 229 891 3 466 018 4 013 231 5 096 922 5 325 358 5 338 247 5 354 623 5 373 151 5 373 325 5 390 008 5 596 326 5 912 456 5 934 873 1 512 986 1 867 038

Start

130 540 450 934 537 620 616 476 875 987 1 785 294 1 822 205 2 229 795 3 465 882 4 013 564 5 097 007 5 325 484 5 338 612 5 354 527 5 373 213 5 373 289 5 390 370 5 596 192 5 912 679 5 934 984 1 512 687 1 867 183

Stop

Range of identified transcriptd

30°C 30°C 30°C 10°C 30°C 30°C 10°C 30°C 30°C 30°C 10°C 30°C 30°C 30°C 10°C 30°C 30°C 30°C 10°C 10°C 30°C 10°C

Condition

sRNA genes had been identified by homology with known sRNA genes in other Pseudomonas species as described by 1 Livny and colleagues (2006); 2 Sonnleitner and colleagues (2008); 3 González and colleagues (2008); 4 Filiatrault and colleagues (2010). * RPKM (reads per kilobase and million; Mortazavi et al., 2008). a. Transcripts of identified sRNA overlaps with transcripts of up-stream ORF. b. Identified sRNA gene is located in a predicted protein coding ORF. c. Coordinates of the predicted sRNAs in the P. putida KT2440 genome are calculated based on sRNA sequences deposited at RFAM (http://rfam.sanger.ac.uk/). d. Coordinates of the predicted sRNAs in the P. putida KT2440 genome are according to identified transcripts based on the Illumina cDNA sequencing results.

4

1

1

2

1

1

4

2

3

1

4

1

1

2

4

4

1

2

sRNA

Ref.

RPKM*

Table 2a. Annotated sRNAs with homology in other Pseudomonas species.

P. putida KT2440 genome update 1311


707 293

869 603

1 028 769

IGR 0752

IGR 0886

Start

IGR 0601

# IGR

1 029 120

870 153

707 527

Stop

351

550

234

Length

Identified cDNA reads

136 135 133 98 212 212 212 261 272 252 272 254

254 254 254 254 254 254 254 254

254 254

254 254 254

36 28 27 27

26 26 26 26 26 26 26 26

26 26

26 26 26

Start

25 (2) 16 16 20 (1) 18 18 18 47

Length

279

279

279

279 279

279

279 279

279 279 279 279 279

298 280

279

307

160 150 148 117 229 229 229 307

Stop

Length/coordinates of the cDNA sequence motif that can hybridize to targeta

5 751 729

5 084 575

4 983 803

4 858 655 4 896 431

4 593 749

3 252 642 4 583 474

91 529 610 797 1 064 991 1 347 032 2 147 142

81 048 5 927 951

1 942 373

1 273 318

3 356 554 2 074 190 4 356 828 4 841 715 2 984 601 4 111 967 4 113 013 4 333 119

Start

5 751 754

5 084 550

4 983 778

4 858 630 4 896 456

4 593 774

3 252 617 4 583 499

91 554 610 822 1 064 966 1 347 057 2 147 167

81 022 5 927 925

1 942 346

1 273 283

3 356 578 2 074 205 4 356 843 4 841 696 2 984 584 4 111 984 4 113 030 4 333 073

Stop

5′ (230 bp), 3′ (67 bp) 5′ (287 bp), 3′ (73 bp) 5′ (176 bp), 5′ (270 bp)

-

+

-

3′ (63 bp) 3′ (65 bp), 5′ (118 bp)

+

PP4391, PP4392 PP4473, PP4474 PP5045, PP5046

PP4274 PP4304, PP4305

PP2850 PP4059, PP4060 PP4067, PP4068

3′d 5′ (50 bp), 3′ (96 bp) 3′ (58 bp), 5′ (216 bp) + +

PP0082 PP0525 PP0922 PP1172 PP1904

3′ (15 bp) Intragenic 3′ (25 bp) Intragenic 3′ (18 bp)

+ + + +

PP0070 PP5195, PP5196

PP2952 PP1850 PP3831 PP4257 PP2611 PP3617 PP3618 PP3804, PP3805 PP1113, PP1114 PP1741

ORF no.

3′ (79 bp) 5′ (40 bp), 3′ (87 bp)

Intragenic Intragenic Intragenic Intragenic Intragenic Intragenic Intragenic 3′ (15 bp), 3′ (23 bp) 3′ (76 bp), 3′ (57 bp) 3′ (27 bp)

Genomic localization of subject sequenceb,c,d,e

-

-

-

+

+, -

+ -

+ + +

+ -

+

+, -

+ + + + +, -

Strand

Target characteristics

-

-

-

+ + + + + -

Strand

Genomic coordinates

Table 2b. Putative sRNAs that carry sequence motifs for potential binding to P. putida KT2440 targets.

Pyridoxal-phosphate dependent enzyme family protein, SEC-C domain-containing protein Substrate-binding region of ABC-type glycine betaine transport system Sua5/YciO/YrdC/YwlC family protein Binding-protein-dependent transport systems inner membrane component, iron ABC transporter, periplasmic iron-binding protein, putative Tryptophan synthase subunit alpha B12 family TonB-dependent receptor Hypothetical Hypothetical UDP-N-acetylenolpyruvoylglucosamine reductase Hypothetical Trehalose synthase, alpha-amylase family protein Acetyl-CoA carboxylase, biotin carboxylase, putative, Cro/CI family transcriptional regulator NAD-dependent DNA ligase LigA Cation transporter, voltage-gated ion channel cation transporter, sulfate ABC transporter, periplasmic sulfate-binding protein Flagellar basal body rod protein FlgB, chemotaxis protein methyltransferase CheR Alanyl-tRNA synthetase, succinylglutamate desuccinylase Thiamine biosynthesis protein ThiI, glutamine synthetase, type I

LysR family transcriptional regulator MFS1 DNA topo isomerase Cbb3-type cytochrome oxidase component Hypothetical Hypothetical Hypothetical Hypothetical

Encoded protein

1312 S. Frank et al.


Start

2 426 959

2 784 791

2 857 768

4 073 623

4 425 079

# IGR

IGR 2128

IGR 2438

IGR2510

IGR 3586

IGR 3917

Table 2b. cont.


4 425 593

4 074 030

2 858 344

2 784 937

2 427 200

Stop

514

407

576

146

241

Length


56 71 72 69 40 45 39 48

353 (1) 338 337 340 (3) 34 (1) 29 35 (2) 26 121 321 45

56 56

353 353 (1)

122 (17) 65 (7) 28 (2)

396 50

43 (5) 359 (1)

8

21 (1) 338 371 396 396

8

21 (1)

159 (2) 58 (4) 56 (6) 33 (2)

203 155 152 2 7

254

26

15 19 (1) 19 (1) 27 22 (1)

Start

Length

242 385 72

73 73

73 73

408 408

408 408

408 408

438 408

496 428 451 428

28

28

217 173 170 28 28

279

Stop


4 301 642 4 301 854 2 674 992

4 678 627 386 715

1 065 043 2 054 689

5 545 316 4 302 318

3 826 241 2 925 774

5 222 562 1 298 543

1 420 153 1 442 064

2 855 970 1 995 842 5 390 520 5 149 120

2 046 732

1 363 011

729 312 3 055 608 3 296 945 3 300 101 1 088 508

5 827 831

Start

4 301 763 4 301 918 2 674 965

4 678 661 386 690

1 065 010 2 054 717

5 545 652 4 302 656

3 826 593 2 925 437

5 222 914 1 298 191

1 420 195 1 441 706

2 856 128 1 995 898 5 390 465 5 149 152

2 046 752

1 362 991

729 298 3 055 590 3 296 927 3 300 127 1 088 529

5 827 806

Stop

Intragenic Intragenic 5′ (63 bp) 3′ (2 bp) 3′ (6 bp) 3′ (36 bp) 3′ (8 bp) 5′ (3 bp) 3′ (138 bp) 5′ (247 bp) 5′ (245 bp), 3′ (308 bp) 5′ (143 bp) 3′ (52 bp), 3′ (22 bp) 5′e 3′ (52 bp), 3′ (25 bp) 3′ (71 bp) 5′ (188 bp), 3′ (40 bp) 3′ (82 bp) 3′ (810 bp), 3′ (496 bp) 3′ (13 bp) 5′ (67 bp), 5′ (140 bp) 3′ (38 bp) 5′ (157 bp), 5′ (144 bp) 3′ (134 bp) 3′ (346 bp) 3′ (219 bp) 3′ (58 bp)

+ + + + + +

+ + -

+ -

+

+

+ -

+ -

+ -

3′ (17 bp), 3′ (157 bp)


-

Strand

Genomic coordinates

PP2508 PP1782 PP4739 PP4534, PP4535 PP1244 PP1260, PP1261 PP4602 PP1133, PP1134 PP3381 PP2563, PP2564 PP4879 PP3774, PP3775 PP0923 PP1829, PP1830 PP4141 PP0321, PP0322 PP3774 PP3774 PP2343, PP2344

PP1819

phoQ

PP0623 PP2667 PP2899 PP2901 PP0942

PP5107, PP5108

ORF no.

+ + +, -

-, +

-, +

+, -

-

+, -

+

+ + -

+

+

+ + + +

+, -

Strand


Phospholipid/glycerol acyltransferase Alpha/beta fold family hydrolase, chorismate synthase DNA polymerase III subunit epsilon Threonine aldolase, serine hydroxymethyltransferase Alginate lyase 2 Alginate lyase 2 Hypothetical

Hypothetical ISPpu9, transposase, 2-hydroxyacid dehydrogenase AraC family transcriptional regulator ISPpu9, transposase, FAD dependent oxidoreductase ISPpu9, transposase Antibiotic biosynthesis protein, putative, Pirin domain protein RNA methyltransferase Alginate lyase 2, sarcosine oxidase, putative

ABC transporter Hypothetical Penicillin amidase family protein pmba protein, putative modulator of DNA gyrase Integral membrane sensor signal transduction histidine kinase Methyl-accepting chemotaxis sensory transducer Hypothetical dTDP-4-dehydrorhamnose 3,5-epimerase Hypothetical Hypothetical

Monofunctional biosynthetic peptidoglycan transglycosylase, RNA polymerase factor sigma-32

Encoded protein


4 630 234

5 047 062

5 148 919

5 391 300

6 016 689

IGR 4095

IGR 4451

IGR 4535

IGR 4740

IGR 5267

6 016 901

5 391 600

5 149 404

5 047 421

4 631 028

4 619 242

Stop

212

300

485

359

794

497

Length


203 11 201 1 202 202 11 15 6 8 132 110 133 44 132

34 (2) 41 (2/1) 33 (2) 33 (2) 30 (2) 34 (3) 35 (4) 20 (1) 26 (2) 30 (1) 25 (2) 20 (1) 26 (3)

47 47 51

20 258

2 238

Start

42 (3) 31 (1)

20 17 16

756 (54) 141 (21)

497 (32) 144 (23)

Length

40 27 157 138 157 63 157

234 40 234 234 40 48

244 41

66 63 66

773 398

498 381

Stop


1 995 019 1 419 278 1 604 496 6 015 789 2 897 748 4 534 680 1 955 234

5 390 521 1 419 943 2 856 028 2 858 163 2 855 831 2 855 150

1 995 867 5 390 716

4 472 924 4 471 691 4 188 022

4 618 746 4 613 518

4 630 253 4 613 516

Start

1 994 985 1 419 259 1 604 521 6 015 818 2 897 724 4 534 661 1 955 209

5 390 488 1 419 983 2 856 060 2 858 195 2 855 860 2 855 117

1 995 908 5 390 686

4 472 943 4 471 707 4 188 037

4 619 501 4 613 658

4 630 749 4 613 659

Stop 3′ (31 bp) Intragenicc 3′d 3′, 5′c 3′ (132 bp) 5′ (47 bp) 3′ (131 bp), 5′ (64 bp) 3′ (128 bp) 5′ (50 bp), 3′ (96 bp) 5′ (246 bp) Intragenic 5′ (71 bp) 3′ (192 bp) 5′ (271 bp) 3′ (114 bp) 5′ (240 bp) Intragenic 3′ (19 bp) 5′ (100 bp) 5′ (49 bp) 5′e 3′ (42 bp)

Strand + + +

+ + +

+ + + + + + -

+ -


Genomic coordinates

PP1781 PP1242 PP1407 PP5266 PP2552 PP4022 PP1752

PP4739 PP1244 PP2508 PP2510 PP2508 PP2507

PP4094 PP4084, PP4085 PP4086 PP4084, PP4085 PP3962 PP3962 PP3683, PP3684 PP1732 PP4739

ORF no.

+ + +

+ + + +

+

+ + +

-

+ +

Strand


Number in brackets indicates mismatches within the sequence motif. Number in brackets indicates the length (bp) between either start or stop codon and the first position of the target sequence. Target sequence overlaps intergenic region. Target sequence overlaps stop codon. Target sequence overlaps start codon.

4 618 745

IGR 4086

a. b. c. d. e.

Start

# IGR

Table 2b. cont.

Hypothetical Hypothetical Hypothetical Hypothetical Hypothetical Short chain dehydrogenase/reductase family oxidoreductase O-acyltransferase, putative Hypothetical Sulfate transporter Acetyl-CoA hydrolase/transferase family protein Aromatic-L-amino-acid decarboxylase Hypothetical Hypothetical

dTDP-4-dehydrorhamnose 3,5-epimerase Hypothetical

Putative lipoprotein Putative lipoprotein Hypothetical

Hypothetical Hypothetical, RHS family protein, putative

Hypothetical Hypothetical, RHS family protein, putative

Encoded protein




Fig. 2. Length distribution of observed 5′ UTRs. The lengths of 5′ UTRs detected at 30°C (A) or 10°C (B) are sorted by rank numbers. The exponential decline indicates a random distribution of 5′ UTR lengths.

The remaining 42 unclassified transcripts (Table S1) also harboured sequence repeats of 14–17 nucleotides in length. However, the evidence that these repeats represent specific targets recognized by the sRNAs is not as strong as for the sRNAs with longer recognition sequences because the 15 mers present in most P. putida KT2440 IGRs are frequent in other parts of the genome. More precisely, a global Blast search using all P. putida IGRs against the KT2440 genome revealed an average of five to seven 15 bp hits per IGR (data not shown). The overrepresentation of 15 mers is a general feature of the P. putida genome and hence a 15 mer repeat in a putative sRNA and further genomic positions does not safely classify as a target for the biological function of the transcript.

Novel ORFs cDNA sequencing revealed 143 yet undescribed ORFs (Table S2), which were all predicted by the in silico gene finders GLIMMER (Salzberg et al., 1998) and/or GeneMark (Lukashin and Borodovsky, 1998). Of these 143 ORFs, 129 ORFs had at least one homologous protein in another species. The remaining 14 ORFs with no known homologue showed either robust expression and/or were differentially regulated at 10°C and 30°C (Table S2), supporting the view that these mRNAs are translated into functional proteins. An encoded function could be ascribed by the criterion of sequence homology to just two ORFs (PP2095.1: homologue of ribosomal modulation factor; PP4350.1: D-alanine–D-alanine ligase). Fourteen ORFs were classified by gene category, but the remaining 127 ORFs were conserved hypotheticals. Thus, the initial annotation of the KT2440 genome overlooked genes of still unknown encoded function.

Validation of operons In silico predicted operons were verified by contigs of Illumina sequence reads that span adjacent ORFs (Table S3a). This kind of analysis was only applicable to those operons that were sufficiently highly expressed at 10°C and/or 30°C in standard minimal medium with succinate as carbon source so that an unequivocal contig of small reads was mapped to the gap between ORFs. The in silico prediction of 154 operons was verified by mapping the cDNA sequence reads onto the KT2440 genome sequence (Table S3a). In contrast, the organization of polycistronic transcripts into operons differed in 57 cases (Table S3b) from that of the most recent version (February 2009) on the Pseudomonas genome web page (http://www.pseudomonas.com) and the DOOR (Database of prOkaryotic OpeRons) web page (Dam et al., 2007; Mao et al., 2009). Remarkable examples are the organization of an ECF sigma factor operon co-transcribed with a siderophore receptor (PP0160PP0163), the inclusion of the elongation factors G and Tu in a ribosomal protein operon (PP0449-PP0452) and the other organization of chemotaxis (PP4338-PP4340) and flagellar biosynthesis operons (PP4374-PP4378, PP4388-PP4391, PP4395-PP4397). Characteristics of 5⬘ untranslated regions (5⬘ UTRs) of mRNA transcripts One hundred and seventy highly expressed single genes or operons that were covered by sequence reads of 300 RPKM or more were suitable for an in silico analysis of the 5′ UTR. The distribution of the lengths of the untranslated region from the transcriptional to the predicted translational start site of mRNAs exponentially declined with length without any preference for a defined length



Fig. 3. Examples of secondary structure prediction for observed 5′ UTRs. The figure shows predicted secondary structures for 5′ UTRs of ORFs PP0741, PP1868 and PP4473. Structure predictions were obtained using the ‘RNAfold WebServer’ provided by the Vienna RNA WebServers (Gruber et al., 2008). Colouring of the bases represents the ‘positioning probability’ of each base at the indicated position, either paired in stem structures or in unpaired positions of the structure.

segment (Fig. 2). This monoexponential decline is compatible with the null hypothesis of a random walk. In other words the length of the 5′ UTR is specific for the individual ORF and no relevant classes of overrepresented 5′ UTR lengths exist in this dataset. The 5′ UTR of eight mRNA transcripts was shorter than 10 nucleotides long. Hence, these transcripts were classified as leaderless (Table S4). Leaderless mRNAs occur in all kingdoms of life, but were assumed to be rather infrequent in Gram-negative bacteria (Moll et al., 2002). Deep cDNA sequencing of the H. pylori transcriptome (Sharma et al., 2010), however, identified 34 leaderless mRNAs among 1576 ORFs (2.1%). Thus the identification of at least eight leaderless mRNAs in P. putida is not an anomalously high number. To search for promoter motifs, the fifty nucleotides upstream of the ATG start were extracted from the genomic sequence and used as input for the web-based motif-finding tool MEME. Surprisingly, neither a TATA box nor the usually conserved -35 bp region were found within this subset of highly expressed genes despite the

fact that both motifs have been reported to occur globally in the related P. aeruginosa genome (Potvin et al., 2008). However, a pentameric A5 was identified in a subset of 44 sites as the singular robust sequence motif. For comparison, the reader may note that the TATA box was detected in 39% of investigated promoters in the H. pylori transcriptome using very similar methodology (Sharma et al., 2010). Next, the 5′ UTR regions were searched with the Web tool ‘RNAfold Server’ (Gruber et al., 2008) for cisregulatory RNA structures in the leader region. The 15 potentially most relevant secondary structures selected by the criteria of thermodynamic stability of individual stem-loop structures and molar proportion in the ensemble of structures are listed in Table S5. The three most stable secondary structures upstream of PP0741, PP1868 and PP4473 are shown in Fig. 3. In summary, about 10% of the examined 5′ UTR contained sequence motifs that can fold into stable stem-loop structures that may function as a cis-regulatory upstream element of gene expression.



Fig. 4. Quantification of expression levels of Illumina cDNA sequencing results. The left diagram shows highly abundant transcripts with normalized read numbers higher than 5000 RPKM (reads per kilobase and million). Genes were defined to be significantly differentially expressed if RPKM10°C > RPKM30°C + 3✓(RPKM30°C) or vice versa. Significantly differentially expressed genes that represent highly abundant transcripts are listed in Table S5. The right diagram displays a zoomed in version with low abundance transcripts expressed less than 10 RPKM. Significantly differentially expressed genes that represent low abundance transcripts with a fold-change ⱖ2.0 are listed in Table S6.

According to our data the length of the 5′ UTR in P. putida is specific for the individual ORF. Neither a preferential length of the 5′ UTR nor any well conserved promoter sequence motifs were identified. Such a genespecific promoter signature that allows a fine-tuning of gene expression by interplay of specific cis-regulatory sequences with the transcription factors in trans may be advantageous for an environmental microorganism like P. putida that is metabolically versatile and can thrive in a broad range of mesophilic habitats. Refined ORFs Deep cDNA sequencing allowed correction of six genes with alternative length, start and stop codons deviating from that predicted in silico. Table S6 compares the initial and the corrected annotation.

favourable feature of transcriptome sequencing is visualized in Fig. 4. The most abundant transcripts are differentiated (left panel) and differences in gene expression at 10°C and 30°C that are small in relative scale, but large in terms of absolute number of sequence reads are clearly discerned. The same argument applies to the least abundant transcripts in the sample (Fig. 4, right panel). The signature of sequence reads at 10°C and 30°C is visible for each individual gene. Tables S7 and S8 list the extremes of the most and least abundant mRNA transcripts. The 50 most abundant mRNA transcripts made up 60% of the total mRNA pool. The most highly expressed genes encode ribosomal proteins, elongation factors, subunits of the ATP synthase or key enzymes of the citric acid cycle (Table S7). In contrast, the categories transporters, porins and outer membrane proteins were overrepresented among the least abundant transcripts (Table S8).

Extreme cDNA copy numbers Deep cDNA sequencing has the advantage that the expression of a gene is inferred from the normalized counts of aligned sequence reads. Sequence reads are individually counted with the same precision irrespective of the total number of reads of the gene in question, whereas hybridization techniques fail to distinguish expression levels within the groups of very lowly and very highly expressed transcripts because of unfavourable background noise or signal saturation respectively. This

Comparison of RNA-seq with RT/PCR kinetics Next, we compared the performance of RNA-seq with that of RT/PCR kinetics to quantify mRNA transcript levels. Two biological replicate pairs of P. putida KT2440 grown at 10°C and 30°C, one pair of which had also been subjected to RNA-seq, were assessed by quantitative RT/PCR kinetics for 11 mRNA transcripts that differed in their absolute abundance and their differential regulation at 10°C and 30°C by up to three orders of magnitude


1318 S. Frank et al. (Fig. S1, Table S9). The primary data demonstrate a slightly larger difference between the RT/PCR comparisons of the two biological replicates (1.1 ⫾ 0.6, mean ⫾ variance) than between the RNA-seq data with the RT/PCR data obtained from the same RNA preparations (1.0 ⫾ 0.6). For both comparisons the mean was close to or matched the expectancy value of 1.0. In other words, the up- and downregulation of genes was detected with similar sensitivity by both techniques. RT/PCR, however, is less sensitive than RNA-seq in the absolute quantitation of transcript. RNA-seq yields data at a linear scale in terms of counts of sequence reads, whereas RT/PCR follows a complex kinetic, exhibits a small window of exponential amplification of template, has a poorer resolution and requires meticulous calibration (Hoof et al., 1991). Comparison of sequencing and microarray transcriptome platforms Aliquots of the same RNA preparation were assessed for their global cDNA expression profile by deep

sequencing (see above) and hybridization onto Affymetrix and Progenika GeneChips. The custom-made Affymetrix high-density oligonucleotide microarray (Ballerstedt et al., 2007; Wierckx et al., 2008) was designed with a pair-wise configuration of 13 perfect match and mismatch 25 mer oligonucleotides per probe set. The probe sets represent 5330 of the 5420 annotated P. putida KT2440 ORFs. Whereas the Affymetrix array allows analysis of single samples, two samples labelled with different dyes are hybridized simultaneously onto the Progenika microarray, which represents each ORF by two adjacent spots of a 50 mer oligonucleotide (Yuste et al., 2006). Figure 5 compares the outcome of transcriptome analysis on the three platforms. ORFs were ranked by signal intensity on the microarrays (Fig. S2) or by normalized sequence reads [number of reads per kb per million reads, RPKM (Mortazavi et al., 2008)]. The Spearman correlation coefficients between the data generated on the Affymetrix and Illumina platforms were 0.74 for the 10°C and 30°C samples (Fig. 5, right panel). RPKM

Fig. 5. Comparison of the three transcriptome platforms using rank number assignment. Rank numbers were assigned to signal intensities or RPKM (reads per kilobase and million) values obtained from the two microarray platforms or the Illumina cDNA sequencing technology, respectively at 30°C and 10°C.



Fig. 6. Rank number evaluation of gene expression. Rank numbers were assigned increasing from highest to lowest expression level. Rank numbers of gene expression for 30°C and 10°C were plotted against one another.

values and processed hybridization data correlated well for the subset of the 2000 most abundant mRNA transcripts, whereas the correlation was rather poor for the remaining 3000 lowly to non-expressed transcripts. Rank numbers established from the Progenika microarray neither correlated with the order of the Affymetrix data nor with the ranking of RPKM values (Fig. 5, left and middle panels). The variable performance of the three technologies probably reflected the differential sensitivity to detect significant differences in gene expression. Figure 6 compares the outcome of transcriptome analysis at the two temperatures. In the case of the Progenika microarrays the rank numbers of individual ORFs showed on average only minor shifts in their position at 10°C and 30°C. The Affymetrix data showed a larger scatter. Extreme shifts in rank numbers were characteristic for the Illumina platform. In other words, the sensitivity to detect differential expression of a gene decreased in the order Illumina > Affymetrix >> Progenika. Of the 159 ORFs that were found to be significantly differentially regulated by all three platforms (Fig. 7), the detected fold-changes were 2.7 (2.2–3.5; 1.7–11.1) [median (inner quartiles; range)] for the Progenika GeneChips, 6.4 (3.9–12.2; 2.1–279) for the Affymetrix Gene Chips and 6.8 (4.1–14.8; 2.1–1231) for the Illumina sequencing platform. Thus, if the least sensitive Progenika array detected a gene to be differentially regulated, the foldchanges were similar for Affymetrix chips and normalized sequence counts. Next, we wanted to know whether this comparable sensitivity of Affymetrix and Illumina technology applied to all genes irrespective of their absolute expression level. For this, the ORFs were sorted by decreasing expression level and the number of differentially expressed genes per

100 comparably expressed ORFs was counted. Illumina sequencing detected differentially expressed genes at the same frequency at all expression levels ranging from the most to the least abundant transcripts (Fig. 8, upper panel). In the case of the Affymetrix microarrays the frequency of significantly differentially expressed genes linearly declined with decreasing absolute expression level (Fig. 8, lower panel). The histogram demonstrates that the sensitivity of the microarray to detect a significant change of expression depended on the absolute signal. The lower the hybridization signal, the less likely the up- or downregulation of a gene could be verified. The ranking of gene expression level correlated between Affymetrix and Illumina platforms – at least for the 30% of the highly expressed genes – but a few sub-

Fig. 7. Venn diagram of significantly expressed genes displaying a differential regulation higher than twofold that were detected with the three transcriptome platforms. Genes were considered to be significantly differentially expressed if calculated P-value and FDR were ⱕ 0.05 and absolute fold-change ⱖ 2.0 for Affymetrix microarrays, if P-value ⱕ 0.05 and FDR ⱕ 0.1 for Progenika microarrays and RPKM10°C > RPKM30°C + 3√(RPKM30°C) (or vice versa) for Illumina cDNA sequencing.



Fig. 8. Intersection of significant regulated genes obtained from either the Affymetrix or the Illumina results into rank number intervals. The numbers of significantly differentially expressed genes detected with either Affymetrix microarrays or Illumina deep RNA sequencing were assigned to Illumina rank number intervals of 100 in steps of 50 (1–100, 51–150, 101–200, etc.). In total 1282 genes were detected with the Affymetrix microarray to be significantly expressed by more than twofold and 2337 genes with the Illumina deep RNA sequencing. Dashed lines indicate the average number of genes per interval (Illumina 32, Affymetrix 23).

stantial exceptions were noted. 126 genes were classified by the GeneChip as not expressed at 10°C and/or 30°C, although a substantial number of cDNA fragments were sequenced by the Genome Analyzer. Table S10 lists

these genes with GeneChip derived rank numbers > 5200 and Illumina derived rank numbers < 3000. Of these 126 genes, 74 genes had not been spotted onto the microarray, 11 of which had been classified as pseudogenes and


P. putida KT2440 genome update 1321 41 as ORFs with a validated frameshift mutation. No sequence reads were found for one further pseudogene and another 17 ORFs containing frameshifts, indicating that no mRNA transcripts exist from these genes under standard growth conditions at 10°C and/or 30°C. Although some genes, such as PP5394, are known to carry authentic frameshift mutations so that the annotation is retained, deep cDNA sequencing provided strong evidence for other coordinates of coding sequence devoid of mutation at eight loci. The updated annotation that is consistent with the experimental cDNA sequence data is listed in Table S11. Of the 126 genes, 52 were represented on the microarray, but no hybridization signal above threshold was detected, although these genes showed a similar distribution of cDNA sequence reads (median RPKM 25; inner quartiles 17–35; range 2.5–358) to the total cDNA sequencing transcriptome (median RPKM 16; inner quartiles 6.4–50; range 0.7–32,415). In other words, the RNAseq. expression profile of the hybridization negative genes is characteristic for P. putida. However, the genes are shorter [mean 514 bp vs. 998 bp (Nelson et al., 2002)] and have lower G + C contents (53.7% vs. 61.6%) than the genomic average (Table S12). It is conceivable that short and GC-poor genes may escape detection under the stringent conditions of microarray hybridization. In other words, an inherent technical bias caused the failure of the microarray to get reliable hybridization signals for this set of short relatively AT-rich genes. Conclusion This study reports on the performance of microarray and cDNA sequencing platforms to assess the transcriptome of P. putida KT2440 at 10°C and 30°C. We would like to draw some conclusions from the obtained data that should globally apply to the field of bacterial genomics and transcriptomics. First, we noted a substantial difference in the sensitivity of the three platforms to uncover significant changes in gene expression. If only a single oligonucleotide per ORF was spotted on the microarray, just a small subset of the differentially regulated genes was detected. The Progenika GeneChip discovered the tip of the iceberg of the differential mRNA expression profiles. The alternative approach, representing each ORF by a probe set of numerous perfect match – single mismatch oligonucleotide pairs, showed a higher sensitivity to detect up- and downregulated genes. cDNA sequencing and Affymetrix microarray data were highly similar within an expression window reaching from the genes with median level to those with high levels of mRNA transcript. However, the sensitivity of the microarray linearly decreased with declining absolute expression levels. In other words, the background noise in terms of

non-specific hybridization prevented the detection of differential expression of many genes that are transcribed at low absolute levels. This bias was not seen in the cDNA sequencing approach. Moreover, the cDNA sequencing technology was superior at displaying shifts in expression levels of the most abundant transcripts. Differential regulation of these top 50 genes escaped detection by microarrays because of signal saturation but was detected with high sensitivity by cDNA sequencing (see Fig. 4, left panel). In summary, cDNA sequencing turned out to be the most unbiased and sensitive technology. cDNA sequencing by Illumina sequencing-by-synthesis technology was not only superior in the detection of up- and downregulated genes than the microarray, but – more importantly – it allowed a re-annotation of the P. putida KT2440 genome. Short cDNA reads were mapped onto the genome sequence. To increase the sensitivity and reliability of the analysis, the transcriptome was studied under two growth conditions that differed in just one easily manageable parameter, the ambient temperature. Thus, the organization of operons and start and stop coordinates of genes were validated or corrected and many novel sRNAs were detected. cDNA sequencing turned out to be a fast and sensitive method of wet-lab validation and correction of the in silico gene predictions. The fast-growing and timely field of small RNAs will particularly benefit from cDNA sequencing because there are still no generally applicable programs for the reliable and comprehensive prediction of sRNAs in bacterial genome sequences. RNA-seq is moreover an efficient method to characterize 5′- and 3′ UTRs at a global level. Whereas the mapping of transcriptional start sites by primer extension protocols is a tedious gene-bygene manoeuvre, cDNA sequencing yields all sites in one experiment assuming that the individual ORF is sufficiently covered by sequence reads. To this end, cDNA sequencing should become an indispensable tool for future projects on re-annotation of sequenced genomes and de novo genome projects of yet unsequenced taxa (Wang et al., 2009). Besides this work on the P. putida genome, transcriptome sequencing has already yielded an improved annotation of the P. syringae (Filiatrault et al., 2010), H. pylori (Sharma et al., 2010), L. monocytogenes (Oliver et al., 2009) and S. solfataricus (Wurtzel et al., 2009) genomes.

Experimental procedures Bacterial strains and growth conditions Pseudomonas putida KT2440 (strain DSM6125) (Bagdasarian et al., 1981) was obtained from the German Resource Centre for Biological Material (DSMZ; Braunschweig, Germany). Bacterial cultures were inoculated from a frozen stock culture of P. putida wild type, and incubated at 30°C for 8 h at 250 r.p.m. in LB medium. An aliquot of 0.2 ml was


1322 S. Frank et al. added to 20 ml M9 medium (Na2HPO4 33.9 g l-1, KH2PO4 15.0 g l-1, NaCl 2.5 g l-1, NH4Cl 5.0 g l-1, MgSO4 2 mM, CaCl2 0.1 mM, FeSO4 ¥ 7 H2O 0.01 mM, pH 6.8) supplemented with 15 mM succinate as sole carbon source in a 100 ml flask and incubated overnight at 30°C. For RNA extraction, bacteria were grown in a 1.5 l batch culture (M9 + 15 mM succinate) using the BioFlo 110 Fermenter (New Brunswick Scientific, Edison, NJ, USA) to ensure constant pH, aeration and agitation. When cultures reached mid-exponential phase (OD600 ~ 0.8) the temperature was decreased from 30°C to 10°C. Three samples, each for parallel RNA extraction, were subsequently taken immediately before temperature downshift (30°C) and 2 h after the media had been cooled to 10°C.

provided by Qiagen, but with self-made buffers for the washing elution steps. Samples were bound to the columns with the commercial binding buffer (Qiagen) and then washed twice with 750 ml washing buffer (potassium phosphate 4 mM, 80% ethanol, pH 8.0) and eluted with 435 ml elution buffer (potassium phosphate, 5 mM, pH 8.5). Labelling efficiency was determined by spectrophotometry at 260, 280, 550 and 650 nm. Equal amounts of cDNA of the samples to be compared on one array were mixed and dried in a vacuum concentrator at 40°C for several hours. For hybridization on Affymetrix microarrays, cDNA samples were fragmented and terminally labelled according to the protocol provided by Affymetrix. The efficiency of the labelling reaction was assessed by a gel-shift assay to prevent hybridization of poorly labelled samples.

RNA isolation Total RNA was extracted according to the protocol provided by Qiagen (RNeasy Mini Kit). Bacterial cultures, three in parallel, were grown until the desired cell density (OD600 ~ 0.8) had been reached. For cell harvest, 2 volumes RNAprotect Bacteria Reagent (Qiagen) were added to one volume bacterial culture and mixed vigorously. The solution was incubated at room temperature for 5 min and immediately centrifuged at 5000 g for 10 min. The supernatant was decanted and if required the bacterial pellet could be stored at -20°C for up to 2 weeks or at -80°C for up to 4 weeks. For cell lysis, the cell pellet was resuspended in a 10% aliquot of the initial sample volume of 1 mg ml-1 lysozyme in 10 mM Tris/HCl, 1 mM EDTA, pH 8.0 and incubated at room temperature for 20 min. Then 1.8 ml RLT buffer (Qiagen) containing 1% b-mercaptoethanol were added and mixed intensively, followed by the addition of 1.2 ml ethanol. The RNA solution was purified using the RNeasy Mini Kit, thereby applying stepwise the total volume to one column. On-column DNase digestion was performed twice for 20 min to ensure the complete removal of genomic DNA. RNA integrity and purity was checked by agarose gel electrophoresis.

cDNA synthesis and labelling cDNA synthesis was performed with a statistically distributed mixture of hexanucleotides as primers (random priming). In this study cDNA was generated from about 10 mg total RNA and pooled from three independent RNA samples using SuperscriptII (Invitrogen) reverse transcriptase according to the manufacturer’s protocols. Purified cDNA samples were split into aliquots for subsequent analysis on Progenika and Affymetrix microarrays and direct high-throughput Illumina sequencing. Samples assessed with the Progenika P. putida genome oligonucleotide microarrays were labelled as described by Yuste and colleagues (2006). The Progenika two-dye array uses the fluorescent dyes cyanine-3 (Cy3) and cyanine-5 (Cy5) (GE Healthcare) whereby both reference and test sample are hybridized onto one array. In this study, 4 arrays per experiment were processed in parallel and each sample was labelled twice with each dye, resulting in eight cDNA samples per experiment. At least 5 mg of total cDNA per sample was used for the labelling reaction. Labelled cDNA samples were purified using the PCR purification columns

Semiquantitative RT/PCR kinetics For PCR kinetics (Bremer et al., 1992), 100 ng of the cDNA was amplified in a 50 ml reaction mixture [5 ml 10¥ reaction buffer (Eurogentec), 3.3 ml 25 mM MgCl2, 1 ml DMSO, 10 ml primer solution (5 mM each), 3 ml dNTPs (2 mM each). 1 U Goldstar-DNA-Polymerase (Eurogentec)]. Aliquots were withdrawn every 2 cycles, separated by electrophoresis and stained with ethidium bromide. The relative amounts Ni and Nj of the template cDNA sequences i and j in the reaction mixture were determined from the titration for the first reaction cycle n when the PCR products became visible by ethidium bromide fluorescence during late exponential phase of PCR according to the equation n ( j) − n (i)

Ni Nj = (1 + R )

,

whereby the efficiency R of the used thermocycler during the exponential phase of PCR was determined to be R = 0.78 ⫾ 0.02 for PCR products of 100–800 bp in length within the interval of reaction cycles 10 < n < 35 (Bremer et al., 1992).

P. putida genomic DNA microarrays, processing and data analysis Array hybridization and processing of primary data were performed at the Array Facility at the Helmholtz Zentrum für Infektionsforschung, Braunschweig. For Progenika genomic DNA microarrays, array slides were incubated in preheated blocking buffer at 42°C and after 45 min washed several times by dipping into ddH2O followed by a final washing step in isopropanol. Slides were dried by centrifugation (2000 ¥ g, 5 min) and subsequently used for hybridization. The desiccated cDNA samples were resuspended in 110 ml hybridization buffer and denatured at 95°C for 2 min. Meanwhile, the array slide was placed in the hybridization chamber and preheated to 42°C. A coverslip was placed on top of the slide with white spacers facing down and the cDNA solution was added slowly in between slide and cover until the array area beneath was completely covered. Samples were hybridized at 42°C for 16 h. Hybridized arrays were treated by several washing steps with subsequent centrifugation (2000 g, 5 min). The cDNA samples prepared for Affymetrix microarrays were hybridized by mixing 55 ml cDNA samples and


P. putida KT2440 genome update 1323 155 ml hybridization mixture with subsequent incubation at 50°C for 16 h at 60 r.p.m. (Affymetrix GeneChip Hybridization Oven 640). Afterwards, the hybridization solution was consecutively replaced by 250 ml buffer A [60 mM sodium phosphate, 0.9 M NaCl, 6 mM EDTA, 0.01 (v/v) Tween 20, pH 7.7] and 250 ml buffer B [MES, 2.6 mM NaCl, 0.01 (v/v) Tween 20]. Thereafter, arrays were transferred to the Affymetrix GeneChip Fluidic Station 450 (program used: Midi_euk-2v3) for subsequent washing and staining sequentially with 600 ml SAPE-, Antibody- and again SAPE-solution. Then, microarrays were scanned using the Affymetrix GeneChip Operation Software Version 1.4 (Affymetrix GeneChip Scanner GCS 30) and then processed with Agilent Scan Control and Feature Extraction Software. Scan images had a resolution of 5 m per pixel. For Progenika and Affymetrix arrays different gene expression features were chosen according to the respective microarray design. For Progenika arrays the two channel gene expression feature was used, for Affymetrix microarrays the one channel feature. Signal intensity per spot was calculated by subtracting background from foreground and provided as median value of all measured pixel per spot. Progenika array data were normalized using the LOWESS (locally weighted linear regression) algorithm to remove intensity-dependent effects that can appear by low signal intensities and to avoid a bias as a result of different label efficiencies of the fluorescent dyes (Quackenbush, 2002). Affymetrix arrays were median scale normalized providing a normalization factor for each array. This factor should be similar between arrays both within an experiment as a control for technical bias or between different experiments to reliably analyse gene expression profiles under different conditions. Raw data were provided by the Array Facility and analysed with the freely available Software BRB ArrayTools (http://linus.nci.nih.gov/pilot/ index.htm). Signal intensities were log2 transformed and median normalized over all arrays of one experiment, and an output file was generated providing mean signal intensities, P-value [two-sample t-test (with random variance model)], FDR (false discovery rate) and fold-change for all represented probe sets. Genes were considered to be significantly differentially expressed if both P ⱕ 0.05 and FDR ⱕ 0.05.

Illumina cDNA sequencing Aliquots of 25 mg cDNA were sequenced with the Genome Analyzer II at GATC Biotech AG (Konstanz). For this, the cDNA was nebulized to generate fragments less than 800 bp long. A terminal ‘A’ was then transferred to the 3′ end and cDNA fragments were ligated to adapters, purified and bridge amplified. 36 cycles of sequencing–by-synthesis were performed for each library with the Genome Analyzer GAII SR generating millions of 36 base reads. Illumina Genome Analyzer Pipeline Version 0.2 software was used to qualify reads passing default signal quality filters (Klockgether et al., 2010). Each position in a sequence read is associated with a quality score Qsolexa reflecting the quality of the respective base call (depending on the measured signal during the sequencing process). For a base X, the relation of its probability value p(X) to the Qsolexa score is given by the formula

Qsolexa = 10 log [p ( X ) 1 − p ( X )]. Therefore, a maximal Qsolexa value of 40 indicates a probability of 99.99% for the called base, or an error probability of 0.01% respectively. Lower Qsolexa values typically appear at the end of reads. If no decisive base call could be made, the respective position in the sequence read was designated as ‘N’ (with Qsolexa set to -0.5 automatically). Sequence reads that passed the default signal quality filter and were not aligned by ELAND (Efficient Large-Scale Alignment of Nucleotide Databases; Illumina, unpublished) to a reference of the P. putida rRNA genes were used for gene expression analysis. The reads were subsequently aligned to the P. putida genome (NC_002947.3) using the Bowtie software package (Langmead et al., 2009). Remaining reads mapped to rRNA were subsequently excluded with a custom PERL script. Four nucleotides were trimmed from the 3′ end of each read and a seed size of 28 bp was used, in which two mismatches were permitted. The quality mismatch sum was 100 and results were output in SAM format (command line: bowtie -t putida -l 28 -e 100 – best – sam -3 4 -n 2 -p 7). A summary table was then generated with the integrative web analysis tool Galaxy (Giardine et al., 2005). The functions ‘coverage’ and ‘join’ were used, respectively, to summarize (i) coverage of each ORF from the P. putida NCBI annotation (version NC_002947.3) in terms of total bp and proportion covered by reads and (ii) the number of reads mapped to each ORF. Reads mapped to ORFs had at least 1 bp overlap with the ORF. The two datasets for 30°C and 10°C differed in the absolute number of both total reads and reads that mapped to the genome. In addition, genes differ greatly in length; therefore, reads were normalized as follows: The ORF length was standardized to 1000 bp and the number of reads to one million reads per experiment (RPKM, see Mortazavi et al., 2008). These analyses were then repeated in the same way for the intergenic regions. Coordinates of intergenic regions were gained using the Galaxy function ‘complement’ on the ORF annotation in Galaxy’s ‘interval’ format (Giardine et al., 2005). Mapped reads were then visualized with the Integrated Genome Browser (Nicol et al., 2009). Genes with RPKM of 5000 or higher and 10 or lower were defined as highly and lowly abundant respectively. Genes were considered to be significantly differentially expressed if RPKM30°C > RPKM10°C + 3√RPKM10°C (or vice versa). The 99% confidence interval for the real value N of a Poissondistributed parameter l is given by N = Nexp ⫾ 3√Nexp, whereby Nexp represents the experimentally determined counts. Full transcriptome data are deposited in accordance with MIAME standards at GEO (http://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc=GSE24176), accession code GSE24176.

Update of annotation The distribution of cDNA reads was compared with the predictions of the initial genome project (Nelson et al., 2002) of the arrangement and size of operons. If inconsistencies between the previous in silico prediction and the experimental sequencing data were noted, the structure of the respec-


1324 S. Frank et al. tive operons was corrected. Sequences of intergenic regions with considerable number of transcripts were analysed with gene prediction software to identify possible protein coding open reading frames. GeneMark (http://exon.biology.gatech. edu/) predictions were performed using GeneMark.hmm 2.4 and the GeneMark.hmm heuristic version, where gene prediction is performed with a pre-computed prediction model based on the reference genome (Lukashin and Borodovsky, 1998) or with a heuristic approach that can be used for gene predictions in novel genomes (Besemer and Borodovsky, 1999) respectively. Glimmer (http://www.cbcb.umd.edu/ software/glimmer/), another gene prediction model, uses interpolated Markov models for the identification of protein coding regions (Salzberg et al., 1998; Delcher et al., 1999). Where genes were predicted, the sequence was searched using Blastx for homologous genes in the NCBI protein database for Bacteria and Archaea (http://blast.ncbi.nlm.nih.gov/ Blast.cgi). Predicted ORFs can vary in length because of the presence of multiple possible start sites. Thus, determination of a putative ORF was based on detecting transcripts that best covered the predicted open reading frame. In the case of small RNAs, sequences of previously verified sRNAs in Pseudomonas species were used to identify transcripts on the basis of sequence similarity. To detect further sRNAs among the novel unclassified transcripts that according to Glimmer or GeneMark do not encode proteins, the presence of sequence repeats elsewhere in the KT2440 genome was queried using Blastn (Altschul et al., 1990). The 5′ UTR regions were extracted from the dataset either manually or with the Galaxy software package (Giardine et al., 2005). Cis-regulatory RNA structures were appraised with the Web tool ‘RNAfold server’ (Gruber et al., 2008) and conserved promoter signals with the tool MEME (Bailey et al., 2006).

Acknowledgements This project was partially supported by the BMBF within the framework of the SysMO consortium, part PSYSMO ‘Towards a quantum increase in performance of P. putida as the cell factory of choice for white biotechnology’, project 3: Key determinants of abiotic stress response of P. putida KT2440. S.F and C.D. were recipients of a predoctoral stipend of the DFG-supported IRTG ‘Pseudomonas: Pathogenicity and Biotechnology’ (GRK 653/3). Microarray analyses were performed using BRB ArrayTools developed by Dr Richard Simon and BRB-Array Tools Development Team.

References Altschul, S.F., Gish, W., Miller, W., Myers, E.W., and Lipman, D.J. (1990) Basic local alignment search tool. J Mol Biol 215: 403–410. Bagdasarian, M., Lurz, R., Rückert, B., Franklin, F.C., Bagdasarian, M.M., Frey, J., et al. (1981) Specificpurpose plasmid cloning vectors. II. Broad host range, high copy number, RSF1010-derived vectors, and a hostvector system for gene cloning in Pseudomonas. Gene 16: 237–247. Bailey, T.L., Williams, N., Misleh, C., and Li, W.W. (2006) MEME: discovering and analyzing DNA and protein

sequence motifs. Nucleic Acids Res 34 (Web Server issue): W369–W373. Ballerstedt, H., Volkers, R.J.M., Mars, A.E., Hallsworth, J.E., Santos, V.A., Puchalka, J., et al. (2007) Genomotyping of Pseudomonas putida strains using P. putida KT2440based high-density DNA microarrays: implications for transcriptomics studies. Appl Microbiol Biotechnol 75: 1133–1142. Besemer, J., and Borodovsky, M. (1999) Heuristic approach to deriving models for gene finding. Nucleic Acids Res 27: 3911–3920. Bremer, S., Hoof, T., Wilke, M., Busche, R., Scholte, B., Riordan, J.R., et al. (1992) Quantitative expression patterns of multidrug-resistance P-glycoprotein (MDR1) and differentially spliced cystic-fibrosis transmembraneconductance regulator mRNA transcripts in human epithelia. Eur J Biochem 206: 137–149. Dam, P., Olman, V., Harris, K., Su, Z., and Xu, Y. (2007) Operon prediction using both genome-specific and general genome information. Nucleic Acids Res 35: 288– 298. Delcher, A.L., Harmon, D., Kasif, S., White, O., and Salzberg, S.L. (1999) Improved microbial gene identification with GLIMMER. Nucleic Acids Res 27: 4636–4641. Dulebohn, D., Choy, J., Sundermeier, T., Okan, N., and Karzai, A.W. (2007) Trans-translation: the tmRNAmediated surveillance mechanism for ribosome rescue, directed protein degradation, and nonstop mRNA decay. Biochemistry 24: 4681–4693. Federal Register (1982) Appendix E, Certified host–vector systems. 47: 17197. Filiatrault, M.J., Stodghill, P.V., Bronstein, P.A., Moll, S., Lindeberg, M., Grills, G., et al. (2010) Transcriptome analysis of Pseudomonas syringae identifies new genes, ncRNAs, and antisense activity. J Bacteriol 192: 2359–2372. Ganesan, H., Rakitianskaia, A.S., Davenport, C.F., Tümmler, B., and Reva, O.N. (2008) The SeqWord Genome Browser: an online tool for the identification and visualization of atypical regions of bacterial genomes through oligonucleotide usage. BMC Bioinformatics 9: 333. Giardine, B., Riemer, C., Hardison, R.C., Burhans, R., Elnitski, L., Shah, P., et al. (2005) Galaxy: a platform for interactive large-scale genome analysis. Genome Res 15: 1451–1455. González, N., Heeb, S., Valverde, C., Kay, E., Reimmann, C., Junier, T., et al. (2008) Genome-wide search reveals a novel GacA-regulated small RNA in Pseudomonas species. BMC Genomics 9: 167. Gooderham, W.J., and Hancock, R.E. (2009) Regulation of virulence and antibiotic resistance by two-component regulatory systems in Pseudomonas aeruginosa. FEMS Microbiol Rev 33: 279–294. Gruber, A.R., Lorenz, R., Bernhart, S.H., Neuböck, R., and Hofacker, I.L. (2008) The Vienna RNA websuite. Nucleic Acids Res 36 (Web Server issue): W70–W74. Hoof, T., Riordan, J.R., and Tümmler, B. (1991) Quantitation of mRNA by the kinetic polymerase chain reaction assay: a tool for monitoring P-glycoprotein gene expression. Anal Biochem 196: 161–169. Klockgether, J., Munder, A., Neugebauer, J., Davenport, C.F., Stanke, F., Larbig, K.D., et al. (2010) Genome


P. putida KT2440 genome update 1325 diversity of Pseudomonas aeruginosa PAO1 laboratory strains. J Bacteriol 192: 1113–1121. Langmead, B., Trapnell, C., Pop, M., and Salzberg, S.L. (2009) Ultrafast and memory efficient alignment of short DNA sequences to the human genome. Genome Biol 10: R25. Lapouge, K., Sineva, E., Lindell, M., Starke, K., Baker, C.S., Babitzke, P., and Haas, D. (2007) Mechanism of hcnA mRNA recognition in the Gac/Rsm signal transduction pathway of Pseudomonas fluorescens. Mol Microbiol 66: 341–356. Livny, J., Brencic, A., Lory, S., and Waldor, M.K. (2006) Identification of 17 Pseudomonas aeruginosa sRNAs and prediction of sRNA-encoding genes in 10 diverse pathogens using the bioinformatic tool sRNAPredict2. Nucleic Acids Res 34: 3484–3493. Lukashin, A.V., and Borodovsky, M. (1998) GeneMark.hmm: new solutions for gene finding. Nucleic Acids Res 26: 1107–1115. Mao, F., Dam, P., Chou, J., Olman, V., and Xu, Y. (2009) DOOR: a database of procaryotic operons. Nucl Acid Res 37: D459–D463. Metzker, M.L. (2010) Sequencing technologies – the next generation. Nat Rev Genet 11: 31–46. Moll, I., Grill, S., Gualerzi, C.O., and Bläsi, U. (2002) Leaderless mRNAs in bacteria: surprises in ribosomal recruitment and translational control. Mol Microbiol 43: 239–246. Morozova, O., Hirst, M., and Marra, M.A. (2009) Applications of new sequencing technologies for transcriptome analysis. Annu Rev Genomics Hum Genet 10: 135–151. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., and Wold, B. (2008) Mapping and quantifying mammalian transcriptomes by CDNA sequencing. Nat Methods 5: 621– 628. Nelson, K.E., Weinel, C., Paulsen, I.T., Dodson, R.J., Hilbert, H., dos Santos, V.A.P.M., et al. (2002) Complete genome sequence and comparative analysis of the metabolically versatile Pseudomonas putida KT2440. Environ Microbiol 4: 799–808. Nicol, J.W., Helt, G.A., Blanchard, S.G., Jr, Raja, A., and Loraine, A.E. (2009) The Integrated Genome Browser: free software for distribution and exploration of genome-scale datasets. Bioinformatic 25: 2730–2731. Oliver, H.F., Orsi, R.H., Ponnala, L., Keich, U., Wang, W., Sun, Q., et al. (2009) Deep RNA sequencing of L. monocytogenes reveals overlapping and extensive stationary phase and sigma B-dependent transcriptomes, including multiple highly transcribed noncoding RNAs. BMC Genomics 10: 641. Potvin, E., Sanschagrin, F., and Levesque, R.C. (2008) Sigma factors in Pseudomonas aeruginosa. FEMS Microbiol Rev 32: 38–55. Quackenbush, J. (2002) Microarray data normalization and transformation. Nat Genet 32 (Suppl.): 496–501. Ramos, J.L., Wasserfallen, A., Rose, K., and Timmis, K.N. (1987) Redesigning metabolic routes: manipulation of TOL plasmid pathway for catabolism of alkylbenzoates. Science 235: 593–596. Regenhardt, D., Heuer, H., Heim, S., Fernandez, D.U., Strömpl, C., Moore, E.R., et al. (2002) Pedigree and taxo-

nomic credentials of Pseudomonas putida strain KT2440. Environ Microbiol 4: 912–915. Salzberg, S.L., Delcher, A.L., Kasif, S., and White, O. (1998) Microbial gene identification using interpolated Markov models. Nucleic Acids Res 26: 544–548. Sharma, C.M., Hoffmann, S., Darfeuille, F., Reignier, J., Findeiss, S., Sittka, A., et al. (2010) The primary transcriptome of the major human pathogen Helicobacter pylori. Nature 11: 250–255. Sonnleitner, E., Sorger-Domenigg, T., Madej, M.J., Findeiss, S., Hackermüller, J., Hüttenhofer, A., et al. (2008) Detection of small RNAs in Pseudomonas aeruginosa by RNomics and structure-based bioinformatic tools. Microbiology 154: 3175–3187. Sonnleitner, E., Abdou, L., and Haas, D. (2009) Small RNA as global regulator of carbon catabolite repression in Pseudomonas aeruginosa. Proc Natl Acad Sci USA 106: 21866–21871. van Vliet, A.H. (2010) Next generation sequencing of microbial transcriptomes: challenges and opportunities. FEMS Microbiol Lett 302: 1–7. Wang, Z., Gerstein, M., and Snyder, M. (2009) RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 10: 57–63. Wierckx, N.J., Ballerstedt, H., de Bont, J.A., de Winde, J.H., Ruijssenaars, H.J., and Wery, J. (2008) Transcriptome analysis of a phenol-producing Pseudomonas putida S12 construct: genetic and physiological basis for improved production. J Bacteriol 190: 2822–2830. Wurtzel, O., Sapra, R., Chen, F., Zhu, Y., Simmons, B.A., and Sorek, R. (2009) A single-base resolution map of an archaeal transcriptome. Genome Res 20: 133–141. Yan, Q., Gao, W., Wu, X.G., and Zhang, L.Q. (2009) Regulation of the PcoI/PcoR quorum-sensing system in Pseudomonas fluorescens 2P24 by the PhoP/PhoQ twocomponent system. Microbiology 155: 124–133. Yoder-Himes, D.R., Chain, P.S., Zhu, Y., Wurtzel, O., Rubin, E.M., Tiedje, J.M., and Sorek, R. (2009) Mapping the Burkholderia cenocepacia niche response via highthroughput sequencing. Proc Natl Acad Sci USA 106: 3976–3981. Yuste, L., Hervás, A.B., Canosa, I., Tobes, R., Jiménez, J.I., Nogales, J., et al. (2006) Growth phase-dependent expression of the Pseudomonas putida KT2440 transcriptional machinery analysed with a genome-wide DNA microarray. Environ Microbiol 8: 165–177.

Supporting information Additional Supporting Information may be found in the online version of this article: Fig. S1. Semiquantitative RT/PCR kinetics of mRNA transcripts in P. putida KT2440 grown at 10°C and 30°C. cDNA was subjected to PCR with aliquots taken at the indicated reaction cycles. The gel-separated PCR products were stained with ethidium bromide (Bremer et al., 1992). Fig. S2. Rank number evaluation of signal intensities for Progenika and Affymetrix microarrays. Rank numbers were assigned to mean and median values obtained from either Progenika or Affymetrix microarrays for the two tested con-


1326 S. Frank et al. ditions. Spearman correlation coefficients indicated a higher scatter of signal intensities detected with the Progenika microarrays. Table S1. Summary of unclassified transcripts encoding either proteins or sRNA genes. Table S2. Expression of newly detected P. putida KT2440 genes at 10°C and 30°C based upon cDNA sequences. Table S3a. Operon organization verified by detected polycistronic transcripts that confirm the current prediction at http:// www.pseudomonas.com Table S3b. Operon organization verified by detected polycistronic transcripts that differ from the current prediction at http://www.pseudomonas.com Table S4. Putative leaderless mRNAs. Table S5. 5′ UTRs with the most stable secondary structures. Table S6. Substitution of open reading frames with incorrect coordinates by experimentally validated cDNA sequences. Table S7. Summary of highly abundant transcripts with an RPKM of 5000 or higher in at least one sample.

Table S8. Summary of low abundance transcripts with less than 10 RPKM in both samples. Table S9. Comparison of quantitation of 11 mRNA transcripts by RNA-seq and semiquantitative RT/PCR kinetics. Table S10. Summary of genes that were not detected with the Affymetrix array but showed sufficient expression with Illumina cDNA sequencing. All these genes were assigned rank numbers ⱖ 5200 for Affymetrix signal intensities and rank numbers ⱕ 3000 for Illumina RPKM values. Table S11. Summary of newly annotated ORFs. The loci had been predicted to be non-protein coding before, but showed sufficient expression in Illumina cDNA sequencing. Table S12. Gene length, GC-content and functional annotation of genes not detected with Affymetrix microarrays. Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.
