To Be or Not to Be a Pseudogene: A Molecular Epidemiological Approach to the mclx Genes and Its Impact in Tuberculosis

RESEARCH ARTICLE To Be or Not to Be a Pseudogene: A Molecular Epidemiological Approach to the mclx Genes and Its Impact in Tuberculosis Catarina Lopes Santos 1,2☯ *, Hanna Nebenzahl-Guimaraes 1,2,3☯ , Marta Vaz Mendes 4 , Dick van Soolingen 3,5 , Margarida Correia-Neves 1,2 1Life and Health Sciences Research Institute (ICVS), School of Health Sciences, University of Minho, Braga, Portugal, 2ICVS/3B’s, PT Government Associate Laboratory, Braga/Guimarães, Portugal, 3National Institute for Public Health and the Environment (RIVM), Bilthoven, Netherlands, 4IBMC—Instituto de Biologia Molecular e Celular, Universidade do Porto, Porto, Portugal, 5Department of Medical Microbiology, Radboud University Nijmegen Medical Centre, Nijmegen, Netherlands ☯These authors contributed equally to this work. *[email protected] Abstract Tuberculosis presents a myriad of symptoms, progression routes and propagation patterns not yet fully understood. Whereas for a long time research has focused solely on the patient immunity and overall susceptibility, it is nowadays widely accepted that the genetic diversity of its causative agent, Mycobacterium tuberculosis, plays a key role in this dynamic. This study focuses on a particular family of genes, the mclxs(Mycobacterium cyclase/LuxR-like genes), which codify for a particular and nearly mycobacterial-exclusive combination of protein domains. mclxs genes were found to be pseudogenized by frameshift-causing insertion (s)/deletion(s) in a considerable number of M.tuberculosis complex strains and clinical isolates. To discern the functional implications of the pseudogenization, we have analysed the pattern of frameshift-causing mutations in a group of M.tuberculosis isolates while taking into account their microbial-, patientand disease-related traits. Our logistic regressionbased analyses have revealed disparate effects associated with the transcriptional inactivation of two mclx genes. In fact, mclx2 (Rv1358) pseudogenization appears to be primarily driven by the microbial phylogenetic background, being mainly related to the Euro-American (EAm) lineage; on the other hand, mclx3 (Rv2488c) presents a higher tendency for pseudogenization among isolates from patients born on the Western Pacific area, and from isolates causing extra-pulmonary infections. These results contribute to the overall knowledge on the biology of M.tuberculosis infection, whereas at the same time launch the necessary basis for the functional assessment of these so far overlooked genes. PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 1/20 OPEN ACCESS Citation: Santos CL, Nebenzahl-Guimaraes H, Mendes MV, van Soolingen D, Correia-Neves M (2015) To Be or Not to Be a Pseudogene: A Molecular Epidemiological Approach to the mclx Genes and Its Impact in Tuberculosis. PLoS ONE 10(6): e0128983. doi:10.1371/journal.pone.0128983 Academic Editor: Igor Mokrousov, St. Petersburg Pasteur Institute, RUSSIAN FEDERATION Received: January 23, 2015 Accepted: May 4, 2015 Published: June 2, 2015 Copyright: © 2015 Santos et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: All relevant data are within the paper and its Supporting Information files. Funding: This work was supported by Fundação para a Ciência e Tecnologia (FCT), Portugal, and cofunded by Programa Operacional Regional do Norte (ON.2—O Novo Norte), Quadro de Referência Estratégico Nacional (QREN), through the Fundo Europeu de Desenvolvimento Regional (FEDER), and from Projeto Estratégico –LA 26 –2013–2014 (PEst-C/SAU/LA0026/2013). H.N.-G. received a personal FCT Grant (SFRH/BD/33902/2209). The funders had no role in study design, data collection Introduction Tuberculosis (TB) is an air-borne contagious disease that remains responsible for high rates of morbidity and mortality worldwide: it is estimated that in 2013, nine million people fell ill with TB and 1.5 million died from it [1]. TB's rate of incidence is declining slowly (1.5% per year in average between 2000 and 2013) [1], which is somehow counter-intuitive given the financial effort put into research and prevention frameworks aiming towards its eradication. The reasons for this halting TB twilight, besides those related with health policies, are the still missing links in the understanding of the disease establishment on its latent or active form and TB transmission. TB has many different facets—from a life-lasting silent infection to an active and potentially deathly disease—and classically most of this variability has been attributed to the hosts’ immune competence. However, and more recently, the role of the etiological agent—Mycobacterium tuberculosis (Mtb)—genotype has been gaining more relevance, as several studies came up demonstrating that small genetic variations in clinical isolates or laboratory strains have a significant impact not only in strict microbial characteristics, such as antibiotic resistance (for instance, [2]), but also in factors related with the host-microorganism relationship dynamics, such as disease progression and/or ability to modulate the host's immune response ([3] and references therein). This work presented here is focused on a family of genes that is almost exclusive of Mycobacterium spp. and particularly abundant in members of the M.tuberculosis complex. Although they have been seldom addressed from a functional point of view, their structure suggests they might be involved in transcriptional regulation and response to quorum sensing (sensu lato) stimuli (i.e., communication within bacterial cells or between bacteria and their hosts) [4]. Its uniqueness relies on codifying for a particular combination of domains: an N-terminal CHD (cyclase homology domain) and a C-terminal LuxR HTH (helix turn helix) domain; depending on the domain identification algorithm, an AAA (ATPases associated with several cellular activities)/NB-ARC domain may also be identified between the other two [4– 6]. Whereas the AAA/NB-ARC domain has the general function of binding and hydrolysing ATP and/or participating in the protein oligomerization ([5] and references therein), the CHD is the catalytic domain of the class III nucleotidyl (adenylyl/guanylyl) cyclases [5,6], and the LuxR is a DNA-binding domain mostly (but not exclusively) known to be associated with quorum sensing transcription modulation [4,7,8]. In the genome of the reference strain Mtb H37Rv, one can find three genes that codify for proteins with this particular domain composition: Rv0386, Rv1358 and Rv2488c [4–6]. For practical reasons, these genes will be referred to as mclx1,mclx2 and mclx3 (from Mycobacterium cyclase/LuxR-like genes), respectively, in the remainder of this manuscript. Taking into account the overall lack of physiological and functional studies on Mclx proteins, they can either be viewed as putative cyclase proteins with additional non-cyclase domain(s) attached [5,6], or as putative transcriptional regulators from the LuxR family that may respond to, or be modulated by, ATP or cAMP [4]. Interestingly, the genome of the reference strain Mtb H37Rv codifies for 16 genes with a CHD domain, some of which are predicted to be transmembranar and others (including those from the Mclx family) predicted to be soluble [5,6]. Both their frequency and their unique diversity in terms of domains composition suggest that the cyclase activity may be a key point in Mtb fitness. In order to complete their function —bind ATP (or, less commonly, GTP) and convert it to the secondary signal cAMP (or cGMP)—nucleotidyl cyclases require a series of conserved residues that have been previously characterized and that are responsible for binding a divalent metal, for stabilizing the transition state species and for selecting and/or attaching the substrate (either ATP or GTP) [5,6]. Interestingly, only the Mclx1 has all the necessary residues for cyclase activity, as Mclx2 lacks one of A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 2/20 and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors have declared that no competing interests exist. the residues necessary for the metal binding and a transition-stabilizing asparagine, Mclx3 lacks a transition-stabilizing asparagine. Moreover, both Mclx2 and Mclx3 seem to lack the substrate selectivity residues [5]. Additionally, the Mclx1 is the only family member that has been functionally characterized: not only was this protein found to have a significant (20%) guanylyl cyclase side activity, besides its adenylyl activity [9], but was also found to have a role in virulence [10]. In fact, Mclx1 was found to be required for a cAMP burst in macrophages upon infection that destabilizes the macrophage immune response: loss of Mclx1 resulted in a reduction in the production of tumor necrosis factor (TNF), and a decrease in the bacterial survival and in the immunopathology in the animal tissues [10]. The general association between the LuxR domain and quorum-sensing mechanisms, together with the lack of some canonical residues for the cyclase activity in the Mclx2 and Mclx3, suggest that these proteins may indeed be quorum-sensing-like transcriptional regulators that respond to or bind to ATP and/or cAMP in an allosteric modulatory fashion. A possible relation between cyclic nucleotides and quorum-sensing (sensu lato) is not new: a dynamic relationship between cyclic nucleotides as signals and quorum-sensing regulatory mechanisms has been observed before, either directly or through CRP (cAMP-binding proteins), in organisms such as Vibrio vulnificus [11], Vibrio fischeri [12], and Vibrio cholerae [13], among others. The study presented here describes an integrative analysis combining genomics and epidemiology and is focused on the clinical consequences of mclx variation. Interestingly, we have found a scattered pattern of pseudogenization among mlcx2 and mclx3 genes, with implications at the level of patients' demographic characteristics and TB clinical manifestations. As so, this report establishes a link between the functionality of the proteins encoded by these two genes and the virulence-related fitness and host adaptation ability of the Mtb. Materials and Methods Screening and alignment of the mclx genes from public available genomes The presence and transcriptional integrity of the mclx genes was analysed in a panel of Mtb complex strains and clinical isolates which genome had been completely sequenced. These organisms were selected from the Genome database of the National Center for Biotechnology Information (NCBI), limiting the search by organism—Mycobacterium tuberculosis complex (taxid: 77643)—and including only genomes with the status "Complete" or "Scaffolds or contigs". This search was performed on the 20 th of March 2014 and yielded 187 organisms. To retrieve the mclx gene sequences from these organisms, each genome or set of scaffolds was used as a reference against which the mclx sequences from the reference strain Mtb H37Rv (Rv0386, Rv1358 and Rv2488c) were mapped, using the assembling tools from Geneious R7.1.4 (Biomatters) [14]. Instances when the putative mclx orthologue spanned more than one scaffold/ contig lead to the elimination of the respective organism from the analyses, as to avoid sequencing misreads potentially attributed to scaffolds/contigs junctions. One hundred and fifty strains/clinical isolates were retained for further analyses: two of Mycobacterium africanum,12 of Mycobacterium bovis, eight of Mycobacterium cannetti, and 128 of Mtb (Table 1). The mclx genes were identified and individually aligned against their reference orthologue from Mtb H37Rv to identify nucleotide substitutions, insertions and deletions, using the ClustalW algorithm [15] available in the Geneious software [14](S1A Fig). For most Mtb strains/clinical isolates, the information on the lineage could be retrieved from the literature and/or information on the genome. For the few cases in which this information could not be found, each genome or set of scaffolds was used as a reference against which the regions of difference (RD) from the reference strain Mtb H37Rv, described by Gagneux A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 3/20 Table 1. Pairwise identity and indels in the mclx genes in a panel of 150 Mtb complex strains/clinical isolates 1 . Organism Lineage mclx 1mclx 2mclx 3 Maf GM041182 WA-2 99.97%—del -1406 (T) 99.90% 99.90% Maf K85 WA-2 99.90%—del -1406 (T); del—2481 C) 99.90% 99.97% Mbv 04–303 - 99.97% 99.90% 100.00% Mbv AF2122/97 - 99.90% 99.90% 100.00% Mbv AN5 - 100.00% 99.90% 100.00% Mbv BCG ATCC 35733 - 100.00% 99.90% 99.97% Mbv BCG ATCC 35740 - 99.90% 99.90% 99.97% Mbv BCG China - 100.00% 99.90% 99.80%—ins—2198–2202 (GGCGG) Mbv BCG Frappier - 100.00% 99.90% 99.97% Mbv BCG Korea 1168P - 100.00% 99.90% 99.97% Mbv BCG Mexico - 100.00% 99.90% 99.97% Mbv BCG Moreau - 100.00% 99.90% 99.97% Mbv BCG Pasteur 1173P2 - 100.00% 99.90% 99.90% Mbv BCG Tokyo 172 - 100.00% 99.90% 99.97% Mcn CIPT 140010059 - 99.90% 99.90% 93.40%—ins—2519–2521 (CCA) Mcn CIPT 140060008 - 99.80% 99.60% 93.50%—ins—2521–2523 (ACC) Mcn CIPT 140070002 - 99.90% 99.90% 98.90% Mcn CIPT 140070005 - 99.70% 99.60% 98.60% Mcn CIPT 140070007 - 99.30% 99.50% 99.00% Mcn CIPT 140070008 - 99.70% 99.90% 99.10% Mcn CIPT 140070013 - 99.30% 98.50% 98.80% Mcn CIPT 140070017 - 99.00% 99.20% 97.00%—ins—2531–2533 (GCC) Mtb '98-R604 INH-RIF-EM' EAm 100.00% 99.90%—ins—840 (T) 99.97%—del—2716 (A) Mtb 02_1987 EAs 99.90% 99.97% 100.00% Mtb 1034 EAs 99.97% 99.97% 100.00% Mtb 210 EAs 99.97% 99.97% 99.97% Mtb 43–16836 IO 100.00% 99.97% 99.90%—ins—1393–1394 (TA) Mtb 7199–99 EAm 100.00% 99.90% 100.00% Mtb BT1 EAs 99.97% 99.97% 99.97% Mtb BT2 EAs 99.97% 99.97% 99.97% Mtb BTB05-552 EAm 100.00% 99.97% 100.00% Mtb BTB05-559 EAm 100.00% 99.97% 100.00% Mtb C EAm 99.50%—ins—3195 (T); ins—3211 (G); ins—3218 (G) 99.90% 99.90%—ins—2919 (G); ins—2992 (G) Mtb CAS/NITR204 EAI 98.90%- 15 del/ 1 ins 99.20%- 13 del 98.70%- 19 del/ 2 ins Mtb CCDC5079 EAs 99.97% 99.90%—del—2523–2525 (CGA) 99.97% Mtb CCDC5180 EAs 99.97% 99.97% 99.97% Mtb CDC1551 EAm 100.00% 99.96% 100.00% Mtb CDC1551A 100.00% 99.97% 100.00% Mtb CTRI-2 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb EAI5 100.00% 99.90% 99.90%—del—1388–1392 (TTGCG) Mtb EAI5/NITR206 IO 99.90% 99.97% 100.00% Mtb EAS054 IO 99.97% 99.97% 99.80%—del—1388–1392 (TTGCG) Mtb F11 EAm 99.97% 99.90%—ins—840 (T) 100.00% Mtb FJ05194 EAs 100.00% 99.90% 100.00% (Continued) A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 4/20 Table 1. (Continued) Organism Lineage mclx 1mclx 2mclx 3 Mtb GuangZ0019 EAm 100.00% 99.97% 100.00% Mtb H37Ra EAm 100.00% 100.00% 100.00% Mtb H37RvCO EAm 100.00% 100.00% 100.00% Mtb HKBS1 EAs 99.97% 99.97% 99.97% Mtb HM EAm*99.97% 99.90%—ins—840 (T) 99.00%—ins—1921–1953 (TG...TG) Mtb HN878 EAs 99.97% 99.97% 99.97% Mtb INS_MDR EAm 100.00% 99.97% 100.00% Mtb INS_SEN EAm 100.00% 99.97% 100.00% Mtb INS_XDR EAm 100.00% 99.97% 100.00% Mtb KZN 1435 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb KZN 4207 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb KZN 605 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb KZN R506 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb KZN V2475 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb MTB-489 ? 99.97% 100.00% 99.97% Mtb NA-A0008 ? 99.97%—del—871 (C) 99.90% 99.80%—del—1388–1392 (TTGCG); del —2044 (C) Mtb NA-A0009 ? 99.90%—del—871 (C); ins—2308– 2310 (AG); del—2901 (C) 99.97% 99.80%—del—1388–1392 (TTGCG); del —2104 (G) Mtb NCGM2209 EAs 99.97% 99.97% 100.00% Mtb OSDD071 EAI 100.00% 99.97% 99.90% Mtb OSDD105 EAm 100.00% 99.90% 99.97%—del—900 (C) Mtb OSDD493 EAm 99.97%—del—352 (C) 99.97% 100.00% Mtb PanR0201 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0202 EAm 100.00% 99.97% 100.00% Mtb PanR0205 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0206 EAm 100.00% 99.90% 100.00% Mtb PanR0207 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0208 EAm 100.00% 99.97% 100.00% Mtb PanR0209 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0304 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0305 EAm 100.00% 99.97% 100.00% Mtb PanR0306 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0307 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0308 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0309 EAm 100.00% 99.97% 100.00% Mtb PanR0313 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0314 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0315 EAm 100.00% 99.90% 100.00% Mtb PanR0316 EAm 100.00% 99.97% 99.97% Mtb PanR0317 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0401 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0402 EAm 100.00% 99.10%—ins—840 (T); del—3466 (G) 99.97% Mtb PanR0403 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0404 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0405 EAm 100.00% 99.90%—ins—840 (T) 99.97% (Continued) A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 5/20 Table 1. (Continued) Organism Lineage mclx 1mclx 2mclx 3 Mtb PanR0409 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0410 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0411 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0412 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0501 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0503 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0505 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0602 EAm 100.00% 96.50%- 5 del/ 4 ins 100.00% Mtb PanR0603 EAm 100.00% 99.90%-ins—840 (T) 99.97% Mtb PanR0604 EAm 100.00% 99.90%-ins—840 (T) 99.97% Mtb PanR0605 EAs 99.97% 99.97% 99.97% Mtb PanR0606 EAs 99.97% 99.97% 99.97% Mtb PanR0607 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0609 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0610 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0611 EAm 99.97% 99.90%—ins—840 (T) 99.97% Mtb PanR0702 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0703 EAm 99.97% 99.90%—ins—840 (T) 99.97% Mtb PanR0704 EAm 99.97% 99.90%—ins—840 (T) 99.97% Mtb PanR0707 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0708 EAm 99.97% 99.90%—ins—840 (T) 99.97% Mtb PanR0801 EAm 100.00% 99.90% 100.00% Mtb PanR0802 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0803 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0804 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0805 EAm 99.97% 99.90%—ins—840 (T) 99.97% Mtb PanR0902 EAm 100.00% 99.90% 100.00% Mtb PanR0903 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PanR0904 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR0906 EAm 99.97% 99.90%—ins—840 (T) 99.97% Mtb PanR0907 EAm 100.00% 99.90% 100.00% Mtb PanR0909 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR1005 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR1006 EAm 100.00% 99.90%—ins—840 (T) 99.97% Mtb PanR1007 EAm 99.97% 99.90%—ins—840 (T) 99.97% Mtb PanR1101 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb PR05 ? 100.00% 99.97% 99.80%—del—1388–1392 (TTGCG) Mtb R1207 EAs 99.90% 99.97% 100.00% Mtb RGTB327 EAm*99.90%—ins—454 (G); ins—684 (A) 99.90%—ins—840 (T); ins—1331 (C) 99.90%—ins—766 (G); ins—2844–2845 (GG) Mtb RGTB423 IO*99.97%—ins—2821 (C) 99.90%—ins—390 (A) 99.70%—ins—1121 (A); ins—2160–2162 (GCC); ins—3017 (G) Mtb S96-129 EAm 100.00% 99.97% 100.00% Mtb Beijing/NITR203 EAs 99.90% 99.90% 99.97% Mtb Erdman = ATCC 35801 EAm 100.00% 99.90% 100.00% Mtb Haarlem EAm 99.10%- 7 ins 99.90% 99.90%—del—2264 (G);del—2362 (G) (Continued) A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 6/20 et al.[16] in the definition of the six phylogeographical lineages, were mapped, and the occurrence of the described long sequence polymorphisms was used for the definition of the lineage (marked with an ""inTable 1). In a few cases the RDs fell in scaffolds/contigs junctions, and therefore the presence of polymorphisms could not be accurately determined, precluding the determination of the lineage (marked with a "?" in Table 1). Screening and alignment of the mclx genes from clinical isolates of an epidemiologically characterized cohort To get an insight into the epidemiological/clinical features that may be associated with the occurrence of the mclx pseudogenization, the mclx transcriptional integrity was analysed in a diversified panel of Mtb clinical isolates. This panel is composed of 140 organisms collected and isolated in the Netherlands from 1993 to 2011, which are fully characterized from an epidemiological and clinical point of view (S1 Table)[17]. Demographic and clinical information, provided by the Registration Committee of the Netherlands Tuberculosis Register (NTR) that approved this retrospective study, were linked to the isolates on the basis of gender, date of birth, year of diagnosis and postal code. For 100 of these clinical isolates, whole-genome sequencing was performed ([17] and Nebenzahl-Guimaraes et al., unpublished data) and the single nucleotide polymorphisms and INDELs (insertions or deletions) were called against the Table 1. (Continued) Organism Lineage mclx 1mclx 2mclx 3 Mtb OSDD515 EAm 100.00% 99.97% 99.97%—del—938 (T) Mtb SUMu001 ? 99.97%—del—952 (G) 100.00% 100.00% Mtb SUMu002 ? 99.97%—del—952 (G) 99.97% 99.97%—ins—2487 (A) Mtb SUMu003 ? 100.00% 99.97% 100.00% Mtb SUMu004 ? 100.00% 99.97% 100.00% Mtb SUMu005 ? 100.00% 99.97% 100.00% Mtb SUMu006 ? 100.00% 99.97% 100.00% Mtb SUMu008 ? 99.97%—del—952 (G) 99.97% 100.00% Mtb SUMu010 ? 100.00% 100.00% 100.00% Mtb SUMu011 ? 100.00% 100.00% 100.00% Mtb SUMu012 ? 99.90%—del—2254 (G); del—2962 (G) 100.00% 100.00% Mtb UM 1072388579 ? 100.00% 99.90% 100.00% Mtb UT205 EAm 100.00% 99.90%—ins—840 (T) 100.00% Mtb W-148 EAs 99.97% 99.90%—del—2523–2525 (CGA) 99.97% Mtb WX3 EAs 99.97% 99.97% 99.97% Mtb X122 EAs 99.97% 99.90%—del—2523–2525 (CGA) 99.97% Mtb XDR1219 EAs 99.97% 99.97% 99.97% Mtb XDR1221 EAs 99.90% 99.97% 100.00% 1 The pairwise identity refers to the % of conserved residues of each gene after aligning it with its orthologue from the reference strain M.tuberculosis H37Rv; the putative pseudogenes are highlighted in bold; indel ocurrences are described by their type (ins, insertion; del, deletion) and by their location (considering an alignment with the reference mclx genes from the Mtb H37Rv strain), except when the total number of isolated occurences exceeds 3, in which case only their frequency is indicated; Maf, Mycobacterium africanum; Mvb, Mycobacterium bovis; Mcn, Mycobacterium cannetti; Mtb, Mycobacterium tuberculosis; WA-2, West-African 2; EAm, Euro-American lineage; EAs, East-Asian lineage; IO, Indo-Oceanic; EAI, East-African Indian lineage *, indicates that the information on the lineage was obtained by aligning the H37Rv RD with the genome/scaffolds of the respective organism. doi:10.1371/journal.pone.0128983.t001 A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 7/20 reference strain H37Rv using Breseq software (version 0.23) with a minimum depth of 15x [18](S1B Fig). SNPs with low-quality evidence (i.e. possible mixed read alignment) or within 5 bp of an INDEL were discarded. The presence of INDELs within the mclx coding regions and their potential to disrupt the open reading frame was evaluated. For the remaining 40 clinical isolates, the mclx genes were PCR-amplified in 3 overlapping fragments, each of 1200 to 1400 base pairs (see S2 Table for information on the primers and PCR conditions), and the purified PCR products were sequenced (by GATC Biotech). Each gene fragment was amplified twice and all fragments were sequenced in both directions. The final sequences were mapped against and aligned with their orthologues from the Mtb H37Rv reference strain, using the Geneious R7.1.4 (Biomatters) [14] software (S1C Fig). All alignments were visually inspected and a conservative approach was applied: whenever the sequencing results failed to converge to an obvious consensus, the gene status (functional/pseudogene) was considered to be unknown. For that reason, a few clinical isolates were not considered in the final analyses (final n= 127). Statistical analyses To envisage the relationship between the clinical and epidemiological features of the TB infection and the status of the mclx genes, a separate logistic regression-based analysis was carried out for mclx2 and mclx3. Firstly, simple binary logistic regressions were performed to identify the significant predictors of each mclx genes status. Afterwards, two multiple binary logistic regressions were performed considering only those considered to be significant predictor variables (p<0.2 in the univariate analysis). Particularly in the case of mclx3, a sequential binary regression model was performed considering three sets of variables (patient-related, microorganism-related and disease-related). All statistical analyses were performed using the IBM SPSS Statistics, version 22 (IBM). Phylogenetic tree construction The FASTA files of publicly available NCBI strains were downloaded and pair aligned with the H37Rv reference sequence (NC_000962.gbk) using MAUVE v2.4.0. Multiple sequence alignments (MSA) using 62 robust SNP markers that have been shown to construct high resolution and reproducible phylogenies [19] were then made for 100 of the clinical isolates and 100 of the publicly available NCBI strains. The MSAs were subsequently used to generate a parsimony-based tree using the DNA parsimony algorithm version 3.69 from the Phylip package. Results A number of mclx genes among the M.tuberculosis complex species are likely pseudogenized The mclx genes were previously identified as a family of genes nearly exclusive to the Mycobacterium genus and particularly abundant among the members of Mtb complex [4]. As a distinctive feature, these genes encode proteins with a particular domain architecture, including a CHD (cyclase homology domain) and a LuxR HTH domain. In order to acquire a better understanding of the distribution of this particular group of genes among the Mtb complex, an alignment-based screening was performed against a diversified panel of organisms with their genome fully-sequenced, including elements from the species M.africanum,M.bovis,M.cannetti and Mtb. Interestingly, this strategy revealed that even though mclx genes have in general an overall high degree of sequence conservation, usually with more than 99% nucleotides identical to their orthologues from the reference strain Mtb H37Rv, in a number of them a few INDELs were present, which caused a disruption in the open reading frame leading to the A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 8/20 accumulation of stop codons and to the truncation of the respective sequence (Table 1). To avoid misinterpretations, and keeping in mind that gene size can vary to a certain extent without necessarily implying loss of function, the following criterion was established: any given mclx gene was considered to be a pseudogene whenever its sequence was truncated in more than 25% of the size of the corresponding Mtb H37Rv reference orthologue. While this criterion would require further functional validation, we consider it appropriate for this initial screening, as it likely minimizes the chance of false positives (i.e., genes that could be considered pseudogenes given a small size variation but that in fact retain their full functionality). Overall, our analysis revealed that the pseudogenization is a rather common phenomenon among the mclx genes: 51.3% of the organisms in the analysed panel have at least one of their mclx genes truncated (Table 1 and Fig 1). As shown in Fig 1, the occurrence of the pseudogenes is not evenly distributed among the three different genes: the number of species with a mclx2 pseudogene is much higher than that of species that suffered pseudogenization in the mclx1 and/or the mclx3. In fact, whereas 7.3% and 8.7% of the organisms in the analysed panel have a pseudogenized mclx1 and mclx3, respectively, 39.3% do so for the mcx2. Accordingly, the percentage of organisms with more than 1 mclx pseudogenized is relatively low (3.3%). Notwithstanding, it should be highlighted that in three (out of the four) organisms with two pseudogenes, the pseudogenization events occur in the mclx1 and mclx3(Table 1 and Fig 1), the two genes with lower pseudogenization occurrence, suggesting that the simultaneous loss of mclx2 and any other mclx may be somehow harmful to the microorganisms. Epidemiological assessment of the mclx pseudogenization In order to evaluate whether there was a relationship between the inactivation of these genes in certain bacterial isolates and the epidemiological and clinical characteristics of the disease caused by those isolates, the presence and functionality of the mclx genes was analysed in a Fig 1. Relative frequency of the organisms harbouring pseudogenes in each of the mclx genes. The percentages are relative to the total number of strains/organisms listed in Table 1, and are stratified according to the three main species: Mtb (Mycobacterium tuberculosis), Mbv (Mycobacterium bovis) and Maf (Mycobacterium africanum). doi:10.1371/journal.pone.0128983.g001 A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 9/20 significant variables buttresses the hypothesis that selective pressures and/or the clinical consequences of inactivating mclx2 or mclx3 are dissimilar. This is in agreement with the previously published functional analysis of the mclx1—Agarwall et al. mutated other adenylyl cyclases besides mclx1, namely mclx3, and for none of them significant effects were noticed in survival upon competing at mouse lungs, nor differences in macrophage cAMP levels compared to the wild-type, as it was for the mclx1, strongly suggesting that they have different functions [10]. The univariate and multivariate binary logistic regression analyses uncovered a number of statistically significant relationships that highlight the potential impact of the mclx genes functionality in aspects related with the microorganism biology and fitness, TB development and patients' demography. In mclx2, the main factor associated with pseudogenization appears to be the microorganism lineage, with all but one organism carrying a pseudogenized gene belonging to the EAm lineage among the panel of strains used for this analysis. Accordingly, in the publicly available genomes there is only one pseudogenization of mclx2 outside the EAm lineage. Age and ethnicity might be significant factors as well, with the multivariate model showing that native Dutch patients present an approximately 6.6-fold higher risk of carrying a pseudogenized mclx2, and younger patients having a decreased OR for these forms. Although this could suggest some degree of adaptation it is also true that 78.6% of the native Dutch isolates actually belong to the EAm lineage, and both age and ethnicity lose significance in the multivariate model controling for lineage. As such, even if age and ethnicity do play a role in the pseudogenization of mclx2, this should be a rather limited one. For mclx3, the results were quite dissimilar from those of mclx2. Significant variables in the multivariate model for mclx3 include patient birth region and disease localization. Given the strong phylogeographical structure of Mtb lineages, the relation between the patients' birth region and mclx3 pseudogenization can either be interpreted as a reflection of a phylogenetic signal (the effect of the microorganisms' lineage by itself cannot be evaluated for this gene, as all pseudogenes are found among strains belonging to the IO lineage in the analysed panel) or due to differences in the individuals from different regions (either genetic or socially/culturallyimplemented). The relationship with disease localization remains consistently significant after correcting for all the other variables with a statistically significant signal in the univariate analysis. This suggests that mclx3 plays a key role in the establishment of a pulmonary infection— and therefore its absence causes an adjustment of the infection to the extra-pulmonary space— or, conversely, that Mclx3 function prevents the infection from spreading. Several risk factors for extra-pulmonary forms of TB have been addressed and characterized previously, both in what comes to the microbial influence [21–23], and also regarding host factors [24]. Concerning genetic microbial features, large INDEL polymorphisms in a phospholipase C gene, plcD, have been significantly associated with extra-pulmonary forms of TB when compared to strains without a plcD interruption [22]. On the other hand, the study of the same kind of mutations occurring in other genes from the same family, plcA,plcB and plcC, failed to show such a correlation [23]. This present study parallels this: mclx3 is strongly and significantly associated with extra-pulmonary TB, but such is not the case for its paralogue mclx2. In what concerns host features, gender, ethnicity and HIV-status have all been found to be significant risk factors for extra-pulmonary forms of TB [24]. Whereas in this study we could not access the patients' ethnicity (concerning race/skin colour), both gender and birth region (a possible proxy) were accounted for in the multivariate model. As a precaution, the data was re-analysed integrating the mclx3 pseudogenization status as a putative risk factor for extra-pulmonary infections (as opposed to strictly pulmonary and/or disseminated infections) and correcting for all host factors previously associated with this form of the disease: age, gender, HIV serological status, birth region and ethnicity. In accordance the pseudogenization of mclx3 appears as an independent and highly significant risk factor for extra-pulmonary TB (S7 Table). Finally, A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 16 / 20 another microbial factor that should be taken into consideration is the microbial lineage. A previous report has demonstrated that, compared to the EAs lineage, the EAm, IO and the EAI lineages are significantly associated with extra-pulmonary forms, even after correcting for relevant host factors [21]. This is particularly important in the context of this study, as almost all mclx3 pseudogenes are found among members of the IO lineage (and actually all of them in the panel used for the regression analysis) and therefore could suggest that the mclx3 pseudogenization is a mere phylogenetic signal. Since the occurrence of mclx3 pseudogenes in the analysed sample is restricted to IO strains, it is not possible to correct for the lineage in the multivariate models. However, previous reports support that the EAI lineage represents a fairly similar risk for extra-pulmonary infections as the IO one [21]. Notwithstanding, in the analysed sample and among the EAI lineage there is only one case (5.6%) of strictly extra-pulmonary infection, a value that deviates significantly from the 13 cases (41.9%) observed for the IO lineage (S8 Table). Although it is not possible to completely disregard the phylogenetic hypothesis, our results suggest that the mclx3 pseudogenization is one of the factors that favor extrapulmonary forms of TB, and its prevalence among the IO lineage justifies that same tendency in this lineage. Conversely, the lack of mclx3 pseudogenes among the studied EAI could help to justify the missing tendency for extra-pulmonary infections in this particular sample. In the context of this manuscript, pseudogene is referred to as any gene whose coding sequence has been abruptly terminated by a large or small INDEL event, leading to the accumulation of stop codons and the precocious ending of the putatively codified peptide. It does not, by any means, reflect a status of overall non-functionality. Pseudogenes can have a number of different functions in the cell. In this context, it is important to highlight the Mycobacterium leprae. Although M.leprae holds a large collection of pseudogenes in its genome (approximately 50%), an interestingly high number of them are actually expressed (43%), and some even vary their expression patterns upon infection or in different leprosy patients, suggesting that they can play a role in the virulence of this microorganism [25–28]. More often than not, this expression occurs from pseudogenes that have stop codons in their reading frames, as is the case of the mclx addressed in this study. Therefore, the mclx pseudogenes should be regarded as potentially functional genes, although codifying smaller proteins/peptides than their orthologues or displaying non-codifying functions. Supporting this hypothesis, one study has previously referred an over-expression of the mclx2 after the induction of an alternative sigma factor (sigF)[29]. Interestingly, the SigF binding site was located within the mclx2 coding region, resulting in a 250 residues-shorter protein. This suggests that shorter versions of at least this mclx may hold an important role under defined conditions. Gene pseudogenization has been often associated with the absence of purifying selective pressures, which allow the accumulation of nucleotide substitutions and INDELs. However, in this case, the high degree of sequence conservation at a nucleotide level suggests otherwise. Frameshift-causing INDELs are sometimes theonlydifferenceinthesequences when compared with their orthologues from the reference genome Mtb H37Rv. Conversely, the mclx genes from the closely related M.cannetti hold a much higher degree of sequence divergence but are not pseudogenized. This suggests that the pseudogenization of the mclx genes is either recent and/or the result of defined selective pressures, as opposed to a longer process of genome erosion in the absence of selection. This work, by describing a family of genes selectively pseudogenized in certain isolates, reinforces the recent trend to complement immunological data with the study of bacterial evolution in order to fully understand—and control—TB. A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 17 / 20 Supporting Information S1 Fig. A. Snapshot of a ClustalW alignment in the Geneious software calling a one basepair insertion in codon 279 of the mclx2 gene in a publicly available NCBI strain. B. Snapshot of the breseq software calling the same insertion in an epidemiologically characterized clinical isolate whose genome has been fully sequenced. Displayed are color-coded Illumina sequencing reads mapping to the H37Rv reference sequence (singled out at the top). C. Snapshot of a ClustalW alignment in the Geneious software calling the same insertion in an epidemiologically characterized clinical isolate in which only the mclx genes have been sequenced. (PDF) S1 Table. Clinical and epidemiological characteristics of the 127 Mtb isolates analysed by binary logistic regression. (PDF) S2 Table. PCR primers and conditions for mclx amplification. (PDF) S3 Table. Multivariate ORs for mclx2 including the "microorganism lineage" as a variable. (PDF) S4 Table. Univariate ORs for mclx2 and mclx3 pseudogenization (clinical isolates with unique RFLP and VNTR). (PDF) S5 Table. Multivariate logistic regression models for mclx2 pseudogenization (clinical isolates with unique RFLP and VNTR). (PDF) S6 Table. Multivariate logistic regression models for mclx3 pseudogenization (clinical isolates with unique RFLP and VNTR). (PDF) S7 Table. Multivariate logistic regression models for extra-pulmonary TB infections. (PDF) S8 Table. Crosstabulation of Mtb lineage vs. local of infection. (PDF) Acknowledgments The authors are grateful to the staff of the RIVM mycobacteriological laboratory for their work on the RFLP and VNTR typing of M.tuberculosis isolates. Author Contributions Conceived and designed the experiments: CLS HN-G MVM MC-N. Performed the experiments: CLS HN-G. Analyzed the data: CLS HN-G DvS MC-N. Contributed reagents/materials/analysis tools: DvS MVM MC-N. Wrote the paper: CLS HN-G. References 1. Global Tuberculosis Report 2014: World Health Organization. 2. Zhang H, Li D, Zhao L, Fleming J, Lin N, Wang T, et al. (2013) Genome sequencing of 161 Mycobacterium tuberculosis isolates from China identifies genes and intergenic regions associated with drug resistance. Nat Genet 45: 1255–1260. doi: 10.1038/ng.2735 PMID: 23995137 A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 18 / 20 3. Coscolla M, Gagneux S (2010) Does M.tuberculosis genomic diversity explain disease diversity? Drug Discov Today Dis Mech 7: e43–e59. PMID: 21076640 4. Santos CL, Correia-Neves M, Moradas-Ferreira P, Mendes MV (2012) A walk into the LuxR regulators of Actinobacteria: phylogenomic distribution and functional diversity. PLoS One 7: e46758. doi: 10. 1371/journal.pone.0046758 PMID: 23056438 5. Shenoy AR, Sivakumar K, Krupa A, Srinivasan N, Visweswariah SS (2004) A survey of nucleotide cyclases in actinobacteria: unique domain organization and expansion of the class III cyclase family in Mycobacterium tuberculosis. Comp Funct Genomics 5: 17–38. doi: 10.1002/cfg.349 PMID: 18629044 6. Shenoy AR, Visweswariah SS (2006) Mycobacterial adenylyl cyclases: Biochemical diversity and structural plasticity. FEBS letters 580: 3344–3352. PMID: 16730005 7. Fuqua C, Winans SC, Greenberg EP (1996) Census and consensus in bacterial ecosystems: the LuxR-LuxI family of quorum-sensing transcriptional regulators. Annu Rev Microbiol 50: 727–751. PMID: 8905097 8. Patankar AV, Gonzalez JE (2009) Orphan LuxR regulators of quorum sensing. FEMS Microbiol Rev 33: 739–756. doi: 10.1111/j.1574-6976.2009.00163.x PMID: 19222586 9. Castro LI, Hermsen C, Schultz JE, Linder JU (2005) Adenylyl cyclase Rv0386 from Mycobacterium tuberculosis H37Rv uses a novel mode for substrate selection. FEBS Journal 272: 3085–3092. PMID: 15955067 10. Agarwal N, Lamichhane G, Gupta R, Nolan S, Bishai WR (2009) Cyclic AMP intoxication of macrophages by a Mycobacterium tuberculosis adenylate cyclase. Nature 460: 98–102. doi: 10.1038/ nature08123 PMID: 19516256 11. Kim SP, Kim CM, Shin SH (2012) Cyclic AMP and cyclic AMP-receptor protein modulate the autoinducer-2-mediated quorum sensing system in Vibrio vulnificus. Curr Microbiol 65: 701–710. doi: 10.1007/ s00284-012-0218-0 PMID: 22961036 12. Lyell NL, Colton DM, Bose JL, Tumen-Velasquez MP, Kimbrough JH, Stabb EV (2013) Cyclic AMP receptor protein regulates pheromone-mediated bioluminescence at multiple levels in Vibrio fischeri ES114. J Bacteriol 195: 5051–5063. doi: 10.1128/JB.00751-13 PMID: 23995643 13. Krasteva PV, Fong JC, Shikuma NJ, Beyhan S, Navarro MV, Yildiz FH, et al. (2010) Vibrio cholerae VpsT regulates matrix production and motility by directly sensing cyclic di-GMP. Science 327: 866– 868. doi: 10.1126/science.1181185 PMID: 20150502 14. Drummond AJ, Ashton B, Buxton S, Cheung M, Cooper A, Duran C, et al. (2011) Geneious v5.4, Available from http://www.geneious.com/. 15. Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, et al. (2007) Clustal W and Clustal X version 2.0. Bioinformatics 23: 2947–2948. PMID: 17846036 16. Gagneux S, DeRiemer K, Van T, Kato-Maeda M, de Jong B, Narayanan S, et al. (2006) Variable hostpathogen compatibility in Mycobacterium tuberculosis. Proc Natl Acad Sci USA 103: 2869–2873. PMID: 16477032 17. Nebenzahl-Guimaraes H, Borgdorff MW, Murray MB, van Soolingen D (2014) A novel approach—the propensity to propagate (PTP) method for controlling for host factors in studying the transmission of Mycobacterium tuberculosis. PLoS One 9: e97816. doi: 10.1371/journal.pone.0097816 PMID: 24849817 18. Barrick JE, Yu DS, Yoon SH, Jeong H, Oh TK, Schneider D, et al. (2009) Genome evolution and adaptation in a long-term experiment with Escherichia coli. Nature 461: 1243–1247. doi: 10.1038/ nature08480 PMID: 19838166 19. Coll F, McNerney R, Guerra-Assuncao JA, Glynn JR, Perdigao J, Viveiros M, et al. (2014) A robust SNP barcode for typing Mycobacterium tuberculosis complex strains. Nature Communcations 5, Article number: 4812. 20. Bentley SD, Comas I, Bryant JM, Walker D, Smith NH, Harris SR, et al. (2012) The genome of Mycobacterium africanum West African 2 reveals a lineage-specific locus and genome erosion common to the M.tuberculosis complex. PLoS Negl Trop Dis 6: e1552. doi: 10.1371/journal.pntd.0001552 PMID: 22389744 21. Click ES, Moonan PK, Winston CA, Cowan LS, Oeltmann JE (2012) Relationship between Mycobacterium tuberculosis phylogenetic lineage and clinical site of tuberculosis. Clin Infect Dis 54: 211–219. doi: 10.1093/cid/cir788 PMID: 22198989 22. Yang Z, Yang D, Kong Y, Zhang L, Marrs CF, Foxman B, et al. (2005) Clinical relevance of Mycobacterium tuberculosis plcD gene mutations. Am J Respir Crit Care Med 171: 1436–1442. PMID: 15805187 23. Kong Y, Cave MD, Yang D, Zhang L, Marrs CF, Foxman B, et al. (2005) Distribution of insertionand deletion-associated genetic polymorphisms among four Mycobacterium tuberculosis phospholipase C genes and associations with extrathoracic tuberculosis: a population-based study. J Clin Microbiol 43: 6048–6053. PMID: 16333097 A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 19 / 20 24. Yang Z, Kong Y, Wilson F, Foxman B, Fowler AH, Marrs CF, et al. (2004) Identification of risk factors for extrapulmonary tuberculosis. Clin Infect Dis 38: 199–205. PMID: 14699451 25. Nakamura K, Akama T, Bang PD, Sekimura S, Tanigawa K, Wu H, et al. (2009) Detection of RNA expression from pseudogenes and non-coding genomic regions of Mycobacterium leprae. Microb Pathog 47: 183–187. doi: 10.1016/j.micpath.2009.06.006 PMID: 19555754 26. Singh P, Cole ST (2011) Mycobacterium leprae: genes, pseudogenes and genetic diversity. Future Microbiol 6: 57–71. doi: 10.2217/fmb.10.153 PMID: 21162636 27. Suzuki K, Nakata N, Bang PD, Ishii N, Makino M (2006) High-level expression of pseudogenes in Mycobacterium leprae. FEMS Microbiol Lett 259: 208–214. PMID: 16734781 28. Williams DL, Slayden RA, Amin A, Martinez AN, Pittman TL, Mira A, et al. (2009) Implications of high level pseudogene transcription in Mycobacterium leprae. BMC Genomics 10: 397. doi: 10.1186/14712164-10-397 PMID: 19706172 29. Hartkoorn RC, Sala C, Uplekar S, Busso P, Rougemont J, Cole ST (2012) Genome-wide Definition of the SigF Regulon in Mycobacterium tuberculosis. J Bacteriol. 194(8):2001–9. doi: 10.1128/JB.0669211 PMID: 22307756 A Molecular Epidemiological Approach to the mclx Genes PLOS ONE | DOI:10.1371/journal.pone.0128983 June 2, 2015 20 / 20