The Genome of a Pathogenic Rhodococcus: Cooptive Virulence Underpinned by Key Gene Acquisitions

1Rhodococcus bacteria belong to the mycolic acid-containing group of actinomycetes together with other major genera such as Corynebacterium, Mycobacterium and Nocardia. The genus Rhodococcus comprises more than 40 species widely distributed in the environment, many with biotechnological applications as diverse as the biodegradation of hydrophobic compounds and xenobiotics, the production of acrylates and bioactive steroids, and fossil fuel desulfurization. The rhodococci also include an animal pathogen, Rhodococcus equi, the genome of which we report here.

2R. equi, a strictly aerobic coccobacillus, is a multihost pathogen that causes purulent infections in various animal species. In horses, it is the etiological agent of “rattles”, a lung disease with a high mortality in foals. R. equi lives in soil, uses manure as growth substrate, and is transmitted by the inhalation of contaminated soil dust or the breath of infected animals. Pathogen ingestion may result in mesenteric lymphadenitis and typhlocolitis, and multiplication in the fecal content of the intestine contributes to dissemination in the environment. R. equi causes chronic pyogranulomatous adenitis in pigs and cattle and severe opportunistic infections in humans, often in HIV-infected and immunosuppressed patients. Human rhodococcal lung infection resembles pulmonary tuberculosis and has a high case-fatality rate,.

3R. equi parasitizes macrophages and, like Mycobacterium tuberculosis (Mtb), replicates within a membrane-bound vacuole. A 80–90 kb virulence plasmid confers the ability to arrest phagosome maturation, survive and proliferate in macrophages in vitro and mouse tissues in vivo, and to cause disease in horses. Virulence-associated protein A (VapA), a major plasmid-encoded surface antigen, is thought to mediate these effects –. The vapA gene is located within a horizontally-acquired pathogenicity island (PAI) together with several other vap genes. Equine, porcine and bovine isolates carry specific virulence plasmid types differing in PAI structure and vap multigene complement, suggesting a role for vap PAI components in R. equi host tropism,.

4Apart from the key role of the plasmid vap PAI, little is known about the pathogenic mechanisms of R. equi. We investigated the biology and virulence of this pathogenic actinomycete by sequencing an analysing the genome of strain 103S, a prototypic clinical isolate. With its dual lifestyle as a soil saprotroph and intracellular parasite, R. equi offers an attractive model for evolutionary genomics studies of niche breadth in Actinobacteria. The comparative genomic analysis of R. equi and closely related environmental rhodococi reported here provides insight into the mechanisms of niche-adaptive genome plasticity and evolution in this bacterial group. The R. equi genome also provides fundamental clues to the shaping of virulence in Actinobacteria.

5The genome of R. equi 103S consists of a circular chromosome of 5,043,170 bp with 4,525 predicted genes (Figure S1) and a circular virulence plasmid of 80,610 bp containing 73 predicted genes. Overall G+C content is 68.76%. Table 1 summarizes the main features of the R. equi genome.

6Orthology analyses (Figure S1) and multiple alignments (Figure 1A) with representative published actinobacterial genomes showed the highest degree of homology and synteny conservation with Rhodococcus jostii RHA1. Next in overall genome similarity was Nocardia farcinica, followed by Mycobacterium spp., whereas Streptomyces coelicolor appeared much more distantly related, consistent with 16S rRNA-derived actinobacterial phylogenies. Some phylogenetic studies have been inconclusive, positioning R. equi either with the nocardiae or rhodococci,. Our genome-wide comparative and phylogenomic analyses indicate this species is a bona fide member of the genus Rhodococcus (Figure 1, Figure S2).

7Interestingly, R. equi has a substantially smaller genome than the soil-restricted versatile biodegrader R. jostii RHA1 (9.7 Mb) and two recently sequenced environmental rhodococci, Rhodococcus erythropolis PR4 (6.89 Mb) and Rhodococcus opacus B4 (8.17 Mb) (see http://www.nite.go.jp/index-e.html). The rhodococcal genomes also differ in structure: R. equi and R. erythropolis have covalently closed chromosomes, whereas those of R. jostii and R. opacus are linear (Table 1, Figure S2). Chromosome topology does not seem to correlate with phylogeny, as R. equi and R. erythropolis belong to different subclades, and the latter is the prototype of the “erythropolis subgroup”, which includes R. opacus. Streptomycetes also have large linear (>8.5 Mb) chromosomes, so linearization appears to have occurred independently in different actinobacterial lineages during evolution, apparently in association with increasing genome size.

8The functional content of the R. equi 103S genome is summarized in Figure S3A. About one quarter of the genome corresponds to coding sequences (CDS) involved in central and intermediate metabolism (n = 1,108) and another quarter corresponds to surface/extracellular proteins (n = 1,073). “Regulators” is the next most populated functional category (n = 464, 10.3%). After adjusting for genome size, the number of membrane-associated proteins is average, but the regulome and secretome are clearly larger than in other Actinobacteria (Figure S4A, S4B, S4C), possibly reflecting specific needs associated with the habitat diversity of R. equi, from soil and feces to the macrophage vacuole. R. equi has 23 two-component regulatory systems, more than twice as many as host-restricted Mtb, and more regulators as a function of genome size than S. coelicolor (Figure S4B). About 29% of the genome encodes products of unknown function. This percentage rises to 44.5% for secreted products (Figure S3B), 13% of which are unique to R.equi. Ortholog comparisons with representative closely related mycolata (R. jostii, N. farcinica and Mtb) showed R. equi to have the highest proportion of species-specific surface/extracellular proteins, consistent with its large secretome. By contrast, R. jostii RHA1 has the largest proportion of unique metabolic genes (Figure S5), consistent with its catabolic versatility. Indeed, R. jostii RHA1 is unique among Actinobacteria in its unusual overrepresentation of metabolic genes (Figure S4D).

9The 5.0 Mb R. equi chromosome contains relatively few pseudogenes (n = 14, Table 1), most associated with horizontally acquired regions (n = 10, including two degenerate DNA mobility genes), consistent with a slow “core” gene decay rate. This suggests that the differences in chromosome size between rhodococci result mainly from genome expansion in environmental species rather than contraction in R. equi.

10We analyzed the paralogous families and local DNA compositional biases to assess the impact of gene duplication (GD) and horizontal gene transfer (HGT) in rhodococcal genome evolution (Tables S1, S2). As expected, both contributed to the chromosome size increase, but with different patterns: linear for GD (i.e. similar percentage of duplicated genes, 32.1, 33.2 and 33.6%, in R. equi, R. erythropolis and R. jostii, respectively), and exponential for HGT (9.5, 14.8 and 19.5%, respectively) (Figure 2). A possible explanation is that HGT involves the simultaneous acquisition of several genes (mean no. of genes per HGT “island” in rhodococci, 8.2 to 10.6). The probability of HGT in rhodococci also increases with chromosome size, as indicated by the mean frequencies of HGT events (1 every 87.0, 67.0 and 54.2 genes in R. equi, R. erythropolis and R. jostii, respectively) (Table S1). Moreover, recently acquired HGT islands, mostly containing “non-adapted” DNA dispensable in the short term in the new host species, are likely to evolve more freely and to tolerate further HGT insertions. This may be the case for two large chromosomal HGT “archipelagos” of ≈90 and 190 Kb in 103S, which probably were generated by an accumulation of HGT events. The mosaic structure of these HGT regions and the diversity of source species, as indicated by reciprocal BLASTP best-hit analysis, suggest that they are a composite of several independent HGT events rather than the result of a single “en-block” acquisition (Figures S1 and S6). Rhodococcal genome expansion also involves a linear increase in the number of paralogous families (with larger numbers of paralogs per family) and non-duplicated genes (Table S2), and an increasing number of unique hypothetical proteins (e.g. 164 in R. equi, 408 in R. jostii). Thus, genome expansion in rhodococci involves greater functional redundancy, diversity and innovation.

11About 20% of R. equi HGT islands (Figure S1) are located close to tRNA genes, suggesting the involvement of phages or integrative plasmids in their acquisition. However, almost no DNA mobilization genes or remnants thereof were found associated with HGT regions, suggesting that the lateral gene acquisitions in the R. equi chromosome are evolutionarily ancient. Most HGT genes (52.5%) probably originated from other Actinobacteria, 3.5% of the best hits were from other bacteria, and 44% had no homologs in the databases. Only four integrase genes, one of them degenerate, and an IS1650-type transposase pseudogene were identified in the 103S chromosome. R. equi seems therefore to be genetically stable in terms of mobile DNA element-mediated rearrangements. DNA mobility genes —mostly associated with HGT regions and increasing in abundance with genome size— are more numerous in environmental rhodococci (Table S3). Thus, increasing genetic flux and plasticity are associated with increasing chromosome size in rhodococci.

12Rhodococcal genome expansion can be largely attributed to extrachromosomal elements. R. equi has a single 80 Kb circular plasmid whereas environmental rhodococci have three to five plasmids, including large linear replicons up to 1,123 Kb in size, accounting for a substantial fraction of the genome (e.g. ≈20% in R. jostii RHA1) (Table 1, Table S3). Thus, as observed for chromosomal HGT DNA, the amount of plasmid increases exponentially with genome size. Indeed, one third of the plasmid DNA was HGT-acquired (32.4%, range 19.35–49.7 vs 14.5%, range 9.5–19.5 for the chromosomes), and plasmids may themselves be considered potentially mobilizable DNA. Rhodococcal plasmids also have a much higher density of DNA mobilization genes (Table S3), pseudogenes (Table 1), unique species-specific genes (mean 44.3±16.0% vs 3.6% to 5.6%), and niche-specific determinants (e.g. the intracellular survival vap PAI in R. equi and 11 of the 26 peripheral aromatic clusters in R. jostii ) than the corresponding chromosomes. Rhodococcal plasmids are therefore clearly under less stringent selection and are key players in rhodococcal genome plasticity and niche adaptability.

13No genes with an obvious role in carbohydrate transport were identified in 103S, consistent with the reported inability of R. equi to utilize sugars, confirmed here by Phenotype MicroArray (PMA) screens and growth experiments in chemically defined mineral medium (MM) (Figure S7A). By contrast, R. jostii, R. erythropolis and R. opacus can grow on carbohydrates – and their genomes encode sugar transporters, including phosphoenolpyruvate-carbohydrate phosphotransferase system (PTS) permeases. Interestingly, the intracellular pathogens R. equi, Mtb and Tropheryma whipplei are the only mesophilic Actinobacteria lacking PTS sugar permeases (Table S4). However, Mtb grows on carbohydrates transported via non-PTS permeases. As the PTS is widespread in Actinobacteria, including nonpathogenic rhodococci and mycobacteria, the absence of PTS components in R. equi, Mtb and the genome-reduced obligate endocellular parasite T. whipplei probably results from gene loss.

14The PMA and MM experiments showed that the only carbon sources used by R. equi 103S were organic acids (acetate, lactate, butyrate, succinate, malate, fumarate; but not pyruvate) and fatty acids (palmitate and the long-chain fatty acid-containing lipids Tween 20, 40 and 80) (Figure S7A). In addition to monocarboxylate and dicarboxylate transporters, the 103S genome encodes an extensive lipid metabolic network, with 36 lipases (16 of which secreted) and many fatty acid β-oxidation enzymes, with 40 acyl-CoA synthetases, 48 putative acyl-CoA dehydrogenases, and 23 enoyl-CoA hydratases/isomerases. Thus, R. equi seems to assimilate carbon principally through lipid metabolism. A mutant in the glyoxylate shunt enzyme isocitrate lyase (REQ38290), required for anaplerosis during growth on fatty acids, has severely impaired intramacrophage replication and virulence, indicating that, as reported for Mtb, lipids are a major growth substrate for R. equi during infection in vivo.

15The 103S genome encodes 21 putative amino acid/oligopeptide transporters, and PMA screens and MM growth assays confirmed that R. equi uses several amino acids (tryptophan, tyrosine, phenylalanine, cysteine, methionine) and dipeptides as sources of nitrogen. However, 103S also has pathways for the de novo synthesis of all essential amino acids, consistent with the ability of R. equi to grow in MM containing only an inorganic nitrogen source (Figure S7A). Thus, R. equi can flexibly adapt to fluctuating conditions of amino-acid availability and grow in amino acid-deficient environments, as typically encountered in the infected host by intracellular pathogens. See Figure 3 for a schematic overview of R. equi 103S nutrition and metabolism.

16R. equi strains cannot grow without thiamine and an analysis of the loci involved in its biosynthesis revealed that thiC is absent from 103S, probably due to an HGT event affecting the thiCD genes (Figure S7B, S7C, S7D). The auxotrophic mutation is probably irrelevant for R. equi in the intestine and manure-rich soil owing to the availability of microbially synthesized thiamine. Host-derived thiamine is also probably available to R. equi during infection.

17We investigated the nutritional and metabolic aspects of rhodococcal niche adaptation by comparing the metabolic network of R. equi with that of R. jostii RHA1, the only other rhodococcal species for which a detailed manually annotated genome is available. RHA1 originated from lindane-contaminated soil and was identified by screening for biodegradative capabilities on multiple aromatic compounds, including polychlorinated biphenyls and steroids. Not surprisingly, its genome has an abundance of aromatic degradation pathways and oxygenases involved in aromatic ring cleavage. R. equi is also soil-dwelling but is primarily isolated from clinical specimens and manure-rich environments, involving clearly different selection criteria and habitat conditions. We used reciprocal best-match BLASTP comparisons to identify the species-specific metabolic gene complements, in which the catabolic specialization is likely concentrated. The related pathogenic Actinobacteria, N. farcinica (which shares a dual soil saprophytic/parasitic lifestyle with R. equi) and Mtb (quasiobligate parasite) were also included in the analyses. R. jostii RHA1 contains a disproportionately larger number of unique metabolic genes than R. equi, N. farcinica and Mtb (n = 1,260 or 47.2% of total metabolic CDS vs only 326 to 375 or 22.9 to 29.2%, respectively) (Figure S8). The oversized metabolic network of RHA1 results from an expansion in the number and gene content of paralogous families (Table S5) and nonparalogous genes (643 CDS in RHA1 vs 209 to 288). Only three of the 29 aromatic gene clusters present in R. jostii were identified in the 103S genome. R. equi therefore has a much smaller metabolic network than, and essentially lacks the vast aromatic catabolome of, R. jostii RHA1.

18R. equi resembles other environmental Actinobacteria in being able to produce oligopeptide secondary metabolites. The 103S genome encodes 11 large non-ribosomal peptide synthetases (NRPS), including three involved in siderophore formation (see below). The only polyketide synthase (REQ02050) is involved in the synthesis of mycolic acids. By contrast, RHA1 has 24 NRPS and seven polyketide synthases. Thus, genome expansion in R. jostii has been accompanied by an extensive amplification of secondary metabolism.

19R. equi reduces nitrates to nitrites through a NarGHIJ nitrate reductase (REQ04200-30). There is also a NirBD nitrite reductase (REQ32900-30), a NarK nitrate/nitrite transporter (REQ32940) and a putative nitric oxide (NO) reductase (REQ03280) (Figure 3). nirBD is conserved in environmental rhodococci whereas narGHIJ and REQ03280 are not, indicating that R. equi is potentially well equipped for anaerobic respiration via denitrification, a useful trait for survival in microaerobic environments, as typically found in necrotic pyogranulomatous tissue, the intestine or manure. A narG mutation has been shown to attenuate R. equi virulence in mice, consistent with the bacteria encountering hypoxic conditions during infection, although this may also reflect defective nitrate assimilation in vivo.

20Intriguingly, R. equi possesses a D-xylulose 5-phosphate (X5P)/D-fructose 6-phoshate (F6P) phosphoketolase (Xfp, REQ21880), the key enzyme of the “Bifidobacterium” F6P shunt, which converts glucose into acetate and pyruvate and is the main hexose fermentation pathway in bifidobacteria. Unexpected fermentative metabolism has been detected in some strictly aerobic bacteria, such as Pseudomonas and Arthrobacter, but no NAD+ (anaerobic)-dependent lactate dehydrogenase or other obvious pyruvate fermentation enzyme was identified in 103S. As R. equi does not use sugars, a catabolic role for the F6P shunt is possible only if fed via gluconeogenesis/glycogenolysis. Alternatively, the F6P shunt may function in reverse (anabolic) mode in R. equi, in parallel to gluconeogenesis, directing excess acetate and glyceraldehyde-3-phosphate (GAP), generated from lipid metabolism, into the pentose phosphate pathway (PPP) (Figure 3). R. equi 103S has a lutABC operon (REQ16290-320), recently implicated in lactate utilization via pyruvate in Bacillus.

21R. equi tolerates a wide pH range, but growth is optimal between pH 8.5 and 10 (Figure S9). This alkaline pH is similar to that of untreated manure, potentially providing a selective advantage for colonization of the farm habitat. The 103S genome encodes a urease (REQ45360-410), an arginine deiminase (REQ11880), an AmiE/F aliphatic amidase/formamidase (REQ26530, next to REQ26520 encoding a UreI-like urea/amide transporter in an HGT island) and other amidases which, by releasing ammonia, may favor R. equi growth in acidic host habitats such as the macrophage vacuole (pH≤5.5), the airways or the intestine (typical pH values in horse, 5.3–5.7 and 6.4–6.7, respectively, ).

22Like other soil bacteria, R. equi encodes a large number of σ factors (21 σ70) and stress proteins (e.g. eight universal stress family proteins [Usp], five cold shock proteins, three heat shock proteins and several Clp proteins). It also synthesizes the ppGpp alarmone involved in adaptation to amino acid starvation. R. equi is transmitted by soil dust in hot, dry weather and must therefore resist low water availability and desiccation-associated oxidative damage. There are two ABC glycine betaine/choline transporters (REQ00540-70 and REQ14620-60), an aquaporin (REQ29580), and genes for the synthesis of an exopolysaccharide (see below) and the osmolytes ectoine (ectABC, REQ07850), hydroxyectoine (ectD, REQ07850) and trehalose (REQ27400-30), potentially important for osmoprotection and water stress tolerance. R. equi is well equipped to face oxidative stress, with four catalases, four superoxide dismutases, six alkyl hydroperoxide reductases and two thiol peroxidases. It also synthesizes the unique actinobacterial redox-storage thiol compound, mycothiol, the antioxidant thioredoxin (REQ47340-50), and the protein-repairing peptide-methionine sulfoxide reductases MsrA (REQ01570) and MsrB (REQ20650). Three homologs of the virulence-associated mycobacterial histone-like protein Lsr2 (one plasmid vap PAI-encoded, REQ03140 and 05980 chromosomal), and a Dps family protein (REQ14900, cotranscribed with REQ14890 encoding a CsbD-like putative stress protein ), may protect against oxidative DNA damage. NO reductase REQ03280 and a putative NO dioxygenase (REQ10890) may confer resistance to nitrosative stress (Figure 3).

23R. equi 103S showed a degree of resistance to many antibiotics in the PMA screens, including 13 aminoglycosides, nine sulfonamides, six tetracyclines, 10 quinolones, 18 β-lactams and chloramphenicol. Standard susceptibility tests confirmed the resistance of 103S to a number of clinically relevant antibiotics (Table S7). This correlates with the presence in 103S of an array of antibiotic resistance determinants, including five aminoglycoside phosphotransferases, 10 β-lactamases and four multidrug efflux systems. Except for β-lactamase REQ26610, none of the resistance genes are associated with HGT regions or DNA mobility genes, suggesting they are ancient traits selected to confer resistance to naturally occurring antimicrobials rather than recent acquisitions associated with the medical use of antibiotics. Soil organisms tend to carry multiple drug resistance determinants, and homologs of most R. equi resistance genes are present in the genomes of environmental rhodococci, at the same chromosomal location in some cases (Figure S10).

24Potential virulence-associated determinants were identified in silico based on (i) homology with known microbial virulence factors, (ii) literature mining for Mtb virulence mechanisms, (iii) automated genome-wide screening for virulence-associated motifs and (iv) systematic inspection of HGT genes, the secretome, and of genes shared with pathogenic actinomycetes but absent from nonpathogenic species.

25The 103S genome harbors three complete mce (mammalian cell entry) clusters. Despite their name, the mechanisms by which these clusters contribute to mycobacterial pathogenesis remain unclear. The mce4 operon from R. jostii and its homolog mce2 in R. equi have recently been shown to mediate cholesterol uptake, consistent with emerging evidence that mce clusters constitute a new subfamily of ABC importers,. The recently reported lack of effect of an mce2 mutation on R. equi survival in cultured macrophages does not exclude a role in cholesterol utilization in vivo or in IFNγ-activated macrophages, as shown for an Mtb mutant in the homologous mce operon. The surface-exposed PE and PPE proteins account for ≈7% of the coding capacity of the Mtb genome due to massive gene duplication, and are thought to play an important role in mycobacterial pathogenesis. The R. equi genome also harbors PE/PPE genes, although only a single copy of each (Figure S11A). They lie adjacent in an operon (REQ01750-60) with the PE gene first, as frequently observed in Mtb, possibly reflecting the functional interdependence of the PE and PPE proteins. REQ35460-550 is identical in structure to ESX-4, one of the five Mtb ESX clusters, and to the single ESX cluster present in Corynebacterium diphtheriae. ESX loci encode two small proteins, ESAT-6 (REQ35460) and CFP-10 (REQ35440), and their type VII secretion apparatus, which also mediates the export of PE and PPE proteins. ESAT-6 and CFP-10 form heterodimeric complexes and are major T-cell antigens and key virulence factors in Mtb. R. equi possesses six mmpL genes, encoding members of the “mycobacterial membrane protein large” family of transmembrane proteins, which are involved in complex lipid and surface-exposed polyketide secretion, cell wall biogenesis and virulence. There are also four Fbp/antigen 85 homologs (REQ01990, 02000, 08890, 20840), involved in Mtb virulence as fibronectin-binding proteins and through their mycolyltransferase activity, required for cord factor formation and integrity of the bacterial envelope.

26A nine-gene HGT island (REQ18350-430) encodes the biogenesis of Flp-subfamily type IVb pili, recently described in Gram-negative bacteria. We confirmed the presence of pilus appendages in 103S (Figure 4). Gene deletion and complementation analysis demonstrated that the identified R. equi pili (Rpl) mediated attachment to macrophages and epithelial cells (P. González et al., manuscript in preparation). The rpl island is absent from environmental rhodococci and is unrelated to the pilus determinants recently identified in Mtb and C. diphtheriae,.

27R. equi is thought to produce capsular material,, and an HGT region encompassing REQ40580-780 contains genes potentially responsible for extracellular polysaccharide synthesis. Two other HGT islands encode sortases, transpeptidases that attach surface proteins covalently to the peptidoglycan and which are important for virulence in Gram-positive bacteria. Both srt islands encode the putative substrates for the sortases (secreted proteins of unknown function) (Figure S11B).

28Several secreted products are putative membrane-damaging or lipid-degrading factors, including a transmembrane protein with a putative hemolysin domain (REQ12980), three cholesterol oxidases (REQ06750, REQ26800, and REQ43910/ChoE ), four “cutinases”/serine esterases (REQ00480, REQ02020, REQ08540, REQ46060) with potential phospholipase A activity, and 16 lipases. REQ34990 encodes a secreted lipoprotein homologous to MBP70 and MPB83, two major mycobacterial antigens strongly expressed in Mycobacterium bovis BCG. The REQ34990 product has a FAS1/BigH3 domain involved in cell adhesion via integrins. There are also homologs of two mycobacterial cytoadhesins, the heparan sulfate-binding hemagglutinin HbhA (Rv0475) involved in Mtb dissemination (REQ38170), and the multifunctional histone-like/laminin- and glycosaminoglycan-binding protein Lbp/Hlp (REQ31340) (Figure 3).

29Iron is essential for microbial growth and the ability to acquire ferric iron from the host is directly related to virulence. Two NRPS, Rbt1/IupS (bimodular, REQ08140-60) and IupU (REQ23810), are involved in the formation of catecholic siderophores or “rhequibactins”. A third NRPS homologous to Mycobacterium smegmatis Fxb (REQ07630) may be involved in the formation of an oligopeptide ferriexochelin-like extracellular siderophore. This “rhequichelin” is probably transported by the iupABC (REQ24080-100)-encoded putative siderophore ABC permease, homologous to the M. smegmatis FxuABC ferriexochelin transporter (Figure 3). The redundancy of iron acquisition systems may explain the lack of effect on virulence of individual iupU, rbt1/iupS and iupABC mutations.

30Only a few species-specific putative virulence loci were found in the 103S genome, all in HGT islands (e.g. the plasmid vap PAI or the chromosomal rpl locus). Most (≈90%) of the potential virulence-related determinants identified in R. equi were present in the environmental Rhodococcus spp. and/or had homologs in nonpathogenic Actinobacteria (Table 2, Table S8). These included orthologs of many experimentally-determined Mtb virulence genes, most of which (≈84%) are conserved among nonpathogenic mycobaceria or have close homologs in environmental actinomycetes (Table S9). The case of the mce, ESX, and PE/PPE loci is illustrative. Initially thought to be Mycobacterium-specific virulence traits, members of these multigene families are present in R. equi and in nonpathogenic rhodococci (Table S8), consistent with growing evidence that they are actually widely distributed among high-G+C gram-positives, whether environmental or pathogenic,,. Notwithstanding that some of the unknown function genes of the 103S genome may encode novel, previously uncharacterized pathogenic traits, these observations are consistent with a scenario in which R. equi virulence largely involves the “appropriation” or cooption of core actinobacterial functions, originally selected in a non-host environment. Gene cooption (also known as preadaptation or exaptation) is a key evolutionary process by which traits that have evolved for one purpose are employed in a new context and acquire new roles, thus allowing rapid adaptive changes –. Cooptive evolution operates through critical modifications in gene expression and function. These changes are particularly feasible in the larger genomes of soil bacteria, with a characteristic profusion of regulators and functionally redundant paralogs,. Without the need for major changes, stress-enduring mechanisms and other housekeeping components, such as the cell envelope mycolic acids or the bacterial metabolic network, may directly contribute to virulence by affording nonspecific resistance or by enabling the organism to feed on host components. We suggest that a few decisive niche (host)-adaptive HGT events in a direct ancestor of R. equi, such as acquisition of the plasmid vap “intramacrophage survival” PAI and the rpl “host colonization” HGT island (Figure 4), triggered the rapid conversion of a “preparasitic” commensal organism into a pathogen via the cooption of preexisting bacterial functions.

31Based on the well-established principle that coexpression with pathogenicity determinants is a strong indicator of involvement in virulence,, we sought to identify novel R. equi virulence-associated chromosomal factors through their coregulation with the plasmid virulence genes. The expression profiles of 103S and an isogenic plasmid-free derivative (103SP−) were compared, using a custom-designed genomic microarray and in vitro conditions known to activate (37°C pH 6.5) or downregulate (30°C pH 8.0) the virulence genes of the plasmid vap PAI,. The plasmid had little effect on the chromosome in vap gene-downregulating conditions, but significantly altered expression was observed for numerous genes in vap gene-activating conditions (n = 88 with ≥2 fold change) (Table S10). Most of the differentially expressed genes (68%) were upregulated in the presence of the plasmid. These data suggest that the virulence plasmid activates the expression of a number of chromosomal genes, but whether this upregulation involves direct, specific (potentially virulence related) interactions or incidental pleiotropic effects is unclear.

32To define the extent and nature of the virulence plasmid-chromosome crosstalk, we subjected the microarray expression data to network analysis. Unlike classical pairwise comparisons, the network approach captures higher-order functional linkages between genes, facilitating the graphic visualization of gene interconnections. It is thus more powerful for biological inference and gene prioritization for experimental validation. Noisy data also tend to be randomly distributed in the network structure. We used BioLayout Express3D, an application that constructs three-dimensional networks from microarray data by measuring the Pearson correlation coefficients between the expression profiles of every gene in the dataset. This is followed by graph clustering using the Markov Clustering (MCL) algorithm to divide the network graph into discrete modules with similar expression profiles. Microarray data of 103S bacteria exposed to various combinations of temperature (20°C, 30°C and 37°C) and pH (5.5, 6.5 and 8) were included in the computations to control for the excessive weight of the variable presence/absence of plasmid and strengthen the correlation analysis.

33Figure 5A shows a network representation of the functional connections detected in the R. equi transcriptome with a Pearson correlation threshold r ≥0.85. The graph model grouped the virulence plasmid genes into two distinct coregulated modules or clusters: one comprised 36 of the 73 plasmid genes, alsmost all from the housekeeping backbone (replication and conjugal transfer functions); the other contained 15 of the 26 vap PAI genes together with a number of chromosomal genes (Table S11). The plasmid housekeeping backbone nodes clustered together outside the main regulation network, reflecting functional independence from the rest of the regulome, as would be expected from the autonomous nature of the extrachromosomal replicon (Figure 5A). This indicates that the graph structure is biologically significant and reflects actual functional relationships, validating the network model. By contrast, the vap PAI nodes were clearly embedded in the network and established multiple connections with chromosomal nodes (Figure 5A, Figure S12A), suggesting that the plasmid virulence genes have undergone a process of regulatory integration with the host R. equi genome. About half of the predicted products of the chromosomal vap PAI-coregulated cluster genes are metabolic enzymes, the others being transcriptional regulators and transporters (Table S11C). Raising the correlation threshold to a highly stringent r≥0.95 disintegrated the network graph into a multitude of discrete, unconnected subgraphs (see Dataset S2). This did not substantially alter the structure of the two plasmid gene-containing clusters, but isolated two chromosomal genes, REQ23860 and REQ23850, as the most significantly and strongly coregulated with the vap PAI genes (Figure 5B), suggesting a direct regulatory interaction.

34The genes from the plasmid backbone cluster were expressed constitutively in the conditions tested, whereas those from the vap PAI-coregulated cluster responded strongly to temperature, with activation at 37°C. Chromosomal genes in this cluster, particularly REQ23860 and REQ23850, displayed the same pattern, with downregulation in 103SP− at 37°C, suggesting that plasmid factors are required for their induction at high temperature (Figure S12B, Table S11B). The vap PAI encodes two transcription factors, VirR (orf4) and an orphan two-component regulator (orf8), both of which have been shown to influence vap gene expression, and could be involved in the observed plasmid-mediated thermoregulation of the vap PAI-coexpressed cluster.

35REQ23860 and REQ23850 null mutants were constructed and tested in J774 macrophages to determine whether the observed coregulation with the plasmid vap PAI correlates with a role in virulence. The plasmidless derivative 103SP−, unable to proliferate intracellularly, was used as an avirulent control. The two mutants had a significantly attenuated capacity to grow in macrophages, restored to wild-type levels upon complementation with the deleted genes (Figure 6), indicating that REQ23860 and REQ23850 are required for optimal intramacrophage proliferation. The mutated genes encode an AroQ (type II) chorismate mutase (CM) and a bifunctional anthranilate synthase (AS) with fused TrpE and TrpG subunits, respectively, two key metabolic enzymes catalyzing the initial committed steps in aromatic amino-acid biosynthesis. CM generates prephenate, the first intermediate in the pathway leading to phenylalanine and tyrosine, whereas AS catalyzes the first reaction in tryptophan biosynthesis. Downstream at the same locus, REQ23840 encodes a prephenate dehydrogenase (Figure 7), which catalyzes the oxidative decarboxylation of prephenate to the tyrosine precursor 4-hydroxyphenylpyruvate. The intracellular growth defect caused by the mutations may therefore be related to a diminished capacity for de novo synthesis of aromatic amino acids. The R. equi genome encodes four other CM enzymes (including one in the vap PAI ) and an additional AS (bipartite, one subunit encoded in a trpECBA operon and the other by a solitary trpG gene elsewhere in the chromosome). Through their coregulation with the plasmid vap PAI, the redundant REQ23850-60-encoded chorismate-utilizing enzymes may be important for R. equi intracellular fitness and full proliferation capacity, by enhancing the de novo supply of aromatic amino acids, which generally appear to be present at limiting concentrations in the in vivo replication niche of intramacrophage vacuole-residing microbial pathogens,,.

36Somewhat counterintuitively for an organism with a dual lifestyle as a soil saprotroph and intracellular parasite, the R. equi genome is significantly smaller than those of environmental rhodococci. This may reflect that the main R. equi habitats –herbivore intestine, manure and animal tissues– provide a richer and more stable environment than the chemically diverse and probably nutrient-scarce environments of the nonpathogenic species. In nutrient-poor conditions, the simultaneous use of all available compounds as sources of carbon and energy may offer a competitive advantage, driving the selection of expanded genomes with greater metabolic versatility,. Indeed, the much larger genome of the polychlorinated biphenyl-biodegrading R. jostii RHA1 encodes a disproportionately large metabolic network, with a wider diversity of paralogous families, unique metabolic genes and catabolic pathways. The relatively small number of pseudogenes and virtual lack of DNA mobilization genes in R. equi suggests that this species has not experienced a sudden evolutionary bottleneck with a concomitant relaxation of selective pressure and increase in mutation fixation. The “coprophilic” and parasitic lifestyle specialization of R. equi seems to result from a “non-traumatic” adaptive process in an organism that, despite having suffered some specific functional losses (e.g. sugar utilization, thiamine synthesis), remains an “average” soil actinomycete with a normal-sized genome under strong selection. The greater genomic complexity of the environmental Rhodococcus spp. may reflect a “multi-substrate” niche specialization necessarily linked to the strict selection criteria —for unusual metabolic versatility— under which these species are generally isolated,. Our analyses show that genome expansion in the environmental rhodococci has involved a linear gain of paralogous genes and an accelerated pattern of gene acquisition through HGT and extrachromosomal replicons, which evolve more rapidly and clearly play a critical role in rhodococcal niche specialization.

37The lipophilic, asaccharolytic metabolic profile and capacity for assimilating inorganic nitrogen may be key traits for proliferation in herbivore intestine and feces, which are rich in volatile fatty acids, and in the macrophage vacuole and chronic pyogranulome, presumably poor in amino acids and rich in membrane-derived lipids,. The potential for anaerobic respiration via denitrification may be critical for survival in the anoxic intestine or, as suggested for Mtb,, in necrotic granulomatous tissue. The inability to use sugars, unique among related actinomycetes, may confer a competitive advantage in the intestine and feces, dominated by carbohydrate-fermenting microbiota generating large amounts of short-chain fatty acids, which R. equi use as main carbon source. Alkalophily is probably an advantage in fresh manure, a major R. equi reservoir. R. equi is also well equipped to survive desiccation, important for dustborne dissemination in hot, dry weather, when rhodococcal foal pneumonia is transmitted,.

38R. equi infections are notoriously difficult to treat due to the intracellular localization of the pathogen, compounded by a lack of susceptibility to antibiotics (e.g. penicillins, cephalosporins, sulfamides, quinolones, tetracyclines, clindamycin, and chloramphenicol) (Table S7 and refs. therein). With its panoply of drug resistance determinants, the 103S genome illustrates how naturally selected resistance traits, typically abundant in soil organisms, may have an important impact on the clinical management of microbial infections.

39Finally, our analyses suggest that the appropriation of preexisting core actinobacterial components and functions are key events in the evolution of rhodococcal virulence. Although the underlying notion may be intuitively apparent when considering, for example, the contribution of housekeeping genes to bacterial virulence, here we are identifying it specifically as “gene cooption”, a key mechanism enabling rapid adaptive evolution and the emergence of new traits –. Underpinned by a few critical “host niche-accessing” HGT events, such as acquisition of the “intracellular survival” plasmid vap PAI or the “cytoadhesion” chromosomal rpl locus, this evolutionary mechanism is likely to have facilitated the rapid conversion of what was probably an animal-associated commensal into the pathogenic R. equi. Given the pervasive distribution of the “virulence-associated” gene pool among nonpathogenic species (Tables S8, S9), the notion of cooptive virulence is possibly applicable to all pathogenic actinomycetes and, indeed, universally to bacterial pathogens. The incorporation of adaptive changes in the regulation of the “appropriated” genes is a key mechanism in genetic cooption. Our genome-wide microarray experiments and transcription network analyses indicate that the plasmid vap PAI, essential for intracellular survival and pathogenicity, has recruited housekeeping genes from the rhodococcal core genome under its regulatory influence. Among these are two chromosomal genes encoding key metabolic enzymes involved in aromatic amino-acid biosynthesis, coexpressed with the virulence genes of the vap PAI in response to an increase in temperature to 37°C (the body temperature of the warm-blooded host). These two metabolic genes are required by R. equi for full proliferation capacity in macrophages, providing supporting experimental evidence for the cooptive nature of R. equi virulence. A cooptive virulence model is consistent with the sporadic isolation of “nonpathogenic” (pre-parasitic) Actinobacteria, including environmental rhodococci (e.g. R. erythropolis ), as causal agents of opportunistic infections. An appreciation of the importance of gene cooption in the acquisition of pathogenicity provides a conceptual framework for better understanding and guiding research into bacterial virulence evolution.

40We sequenced the original stock of the foal clinical isolate 103, designated clone 103S, to avoid mutations associated with prolonged subculturing in vitro. Strain 103 belongs to one of the two major R. equi genogroups (DNA macrorestriction analysis, unpublished data), is genetically manipulable, and is regularly used for virulence studies,. Random genomic libraries in pUC19 were pair-end sequenced using dye terminator chemistry on ABI3700 instruments, with subsequent manual gap closure of shotgun assemblies and sequence finishing, as previously described. The 103S genome sequence was manually curated and annotated with the software and databases listed in Table S12. A conservative annotation approach was used to limit informational noise. For phylogenomic analyses, putative core ortholog genes were identified by reciprocal FASTA using a minimum cutoff of 50% amino acid similarity over 80% or more of the sequence. A similarity distance matrix was built with the average percentage amino acid sequence identity obtained by pairwise BLASTP comparisons (distance = 100 − average percent identity of 665 loci) and used to infer a neighbor-joining tree with the Phylip package. The accession numbers of the genome sequences used in comparative analyses are listed in Table S13.

41The sequence from the R. equi 103S genome has been deposited in the EMBL/GenBank database under accession no. FN563149.

42The nutritional and metabolic profile of R. equi 103S and its susceptibility to various drugs were analysed in Phenotype MicroArray screens (Biolog Inc., http://www.biolog.com). Substrate utilization was validated in supplemented mineral medium (MM) containing salts, trace elements, and ammonium chloride as the sole nitrogen source (see Figure S7). For electron microscopy, a bacterial cell suspension in 0.1 M Tris-HCl (pH 7.5) was negatively stained with 1% uranyl acetate and observed at 80.0 kV in a Phillips CM120 BioTwin instrument (University of Edinburgh). Fluorescence microscopy was carried out on paraformaldehyde-fixed bacteria with an R. equi whole-cell rabbit polyclonal antiserum and Alexa Fluor 488-conjugated secondary antibodies (both diluted 1∶1000 in 0.1% BSA).

43Total RNA was obtained from logarithmically growing R. equi bacteria (OD600 = 0.8) in Luria-Bertani (LB) medium, by homogenization in guanidinium thiocyanate-phenol-chloroform (Tri reagent, Sigma) with FastPrep-24 lysing matrix and a FastPrep apparatus (MP bio), followed by chloroform-isopropanol extraction, DNAase treatment (Turbo DNA-free, Ambion) and purification with RNeasy kit (Qiagen). RNA quantity and quality were determined with a Nanodrop (Thermo Scientific) and 2100 Bioanalyzer with RNA 6000 Nano assay (Agilent). RNA samples (500 ng) were amplified with the MessageAmp II-bacteria kit and 5-(3-amionallyl)-UTP (Ambion), labeled with Cy3 or Cy5 NHS-ester reactive dyes (GE Healthcare), and purified with RNeasy MinElute (Qiagen). Whole-genome 8×15K custom microarrays with up to four different 60-mer oligonucleotides per CDS (13,823 probes for the chromosome, 201 for the virulence plasmid) (Agilent) were hybridized in Surehyb DNA chambers (Agilent) with 300 ng of Cy3/Cy5-labeled aRNA, using Gene Expression Hybridisation and Wash Buffer kits (Agilent). Three experimental replicates per condition were analyzed, one with dye swap. The hybridization signals were captured and linear intensity-normalized, with Agilent's DNA microarray scanner and Feature Extraction software. Data were subsequently LOESS-normalized by intensity and probe location and analyzed with Genespring GX 10 software (Agilent). Network analysis of microarray expression data was carried out with Biolayout Express3D 3.0 software, using log base 2 normalized ratios of Cy3/Cy5 signals and methods described in detail elsewhere. Biolayout Express3D is freely available at http://www.biolayout.org/.

44In-frame deletion mutants of REQ23860 and REQ23850 were constructed by homologous recombination, using the suicide vector pSelAct for positive selection of double recombinants on 5-fluorocytosine (5-FC). Briefly, oligonucleotide primer pairs CMDEL1/CMDEL2 and CMDEL3/CMDEL4 were used for PCR amplification of two DNA fragments of ≈1.5 Kb corresponding to the seven 3′- and six 5′-terminal codons plus adjacent downstream and upstream regions of REQ23860. The CMDEL2 and CMDEL3 primers are complementary and were used to join the two amplicons by overlap extension. The PCR product carrying the ΔREQ23860 allele was inserted into pSelAct, using SpeI and XbaI restriction sites; the resulting plasmid was introduced into 103S by electroporation and transformants were selected on LB agar supplemented with 80 µg/ml apramycin. The same procedure was followed for ΔREQ23860, with primers ASDEL 1 to 4. Allelic exchange double recombinants were selected as previously described,. For complementation, the REQ23860-50 genes plus the entire upstream intergenic region were amplified by PCR with CACOMP1 and 2 primers and stably inserted into the R. equi chromosome, using the integrative vector pSET152. PCR was carried out with high-fidelity PfuUltra II fusion HS DNA polymerase (Stratagene). The primers used are shown in Table S14.

45Low-passage (<20) J774A.1 macrophages (ATCC) were cultured in 24-well plates at 37°C, under 5% CO2 atmosphere, in DMEM supplemented with 2mM L-glutamine (Gibco) and 10% fetal bovine serum (Lonza) until confluence (≈2×105 cells/well). J774A.1 monolayers were inoculated at 10∶1 MOI with washed R. equi from an exponential culture at 37°C in brain-heart infusion (BHI, OD600≈1.0). Infected cell monolayers were immediately centrifuged for 3 min at 172×g and room temperature, incubated for 45 min at 37°C, washed three times with Dulbecco's PBS to remove nonadherent bacteria, and incubated in DMEM supplemented with 5µg/µl vancomycin to prevent extracellular growth. After 1 h of incubation with vancomycin (t = 0) and at specified time points thereafter, cell monolayers were washed twice with PBS, detached with a rubber policeman and lysed by incutation for 3 min with 0.1% Triton X-100. Intracellular bacterial counts were determined by plating appropriate dilutions of cell lysates onto BHI. The presence of the virulence plasmid was checked by PCR on a random selection of colonies, using traA- and vapA- specific primers to exclude the possibility of intracellular growth defects being due to plasmid loss. As the intracellular bacterial population at a given time point depends on initial numbers, bacterial intracellular kinetics data are expressed as a normalized “Intracellular Growth Coefficient” according to the formula IGC = (IBt = n−IBt = 0)/IBt = 0, where IBt = n and IBt = 0 are the intracellular bacterial numbers at a specific time point, t = n, and t = 0, respectively.