The Pan-Genome of the Zoonotic Neglected Pathogen Bartonella henselae Reveals Two Groups with Different Patterns of Adaptation to Hosts

Abstract

The genus Bartonella is comprised of Gram-negative re-emerging bacteria like Bartonella henselae, which mainly infects humans and survives inside erythrocytes. This species is transmitted by scratches and bites from domestic cats and usually causes a symptomatic infection in humans, known as Cat-Scratch Disease (CSD). The disease causes multiple clinical signs in humans, such as dermatic, cardiovascular, lymphatic, hepatic and nervous system diseases in immunosuppressed individuals. Although the bacteria are highly relevant for its zoonotic importance worldwide, few studies aimed at characterizing these species genomes and there is still no pan-genome study available. Here, we performed phylogenomic, pan-genome and genome plasticity analyses to determine the epidemiological aspects, the size of the pangenome and its variability in the identified pathogenicity and resistance islands. Altogether, our results showed that the genomes are highly similar, with an almost closed pan-genome. Also, we found two subsets of genomes, composed of 7 and 17 genomes of bacteria. Our results point to the need of sequencing more genomes worldwide to better characterize these variations in the pangenome and understand the patterns of adaptation of this species. The highly conserved genomes from this species are very important for the development of new vaccines and analyses of drug targets against this pathogen. Furthermore, these data may then be used in future works, which will be highly relevant for containing the disease worldwide.

Keywords: Bartonella henselae; Cat-Scratch Disease; Pan-Genome; Genome Plasticity; Phylogenomics

Abbreviations

ACT (Artemis Comparison Tool); B. apis (Bartonella apis); B. henselae (Bartonella henselae); BRIG (BLAST Ring Image Generator); CDSs (Coding DNA Sequences); COG (Cluster of Orthologous Groups); CSD (Cat-Scratch Disease); DNA (deoxyribonucleic acid); GC-content (guanine-cytosine content); GIPSy (Genomic Island Prediction Software); LCBs (Locally Collinear Blocks); Mb (mega base pairs); kb (kilo base pairs); MCL (Markov Clustering); NCBI (National Center for Biotechnology Information); PAIs (Pathogenicity Islands); RIs (Resistance Islands); T4SS (type IV secretion system); USA (United States of America); VFDB (Virulence Factor Database); wgMLST (whole genome Multilocus Sequence Typing).

Introduction

The genus Bartonella is comprised of Gram-negative, fastidious, intracellular and reemerging bacteria [1]. Some species from this genus mainly infect cats, but is also an opportunistic pathogen for humans and are frequently acquired through hematophagous arthropod vectors. Those species are able to infect and survive inside erythrocytes, through a long intra-erythrocytic and intra-endothelial infection, which results in recidivist bacteremia both in humans and other mammals [2]. The genus presents at least 13 human pathogenic species [3], among which the 3 most relevant ones being: Bartonella bacilliformis, which causes the Oroya fever [4]; Bartonella quintana, causing trench fever [5]; and, Bartonella henselae, the causative agent of the Cat-Scratch Disease (CSD).

Domestic cats are the main reservoirs of B. henselae and both immunocompetent and immunosuppressed cats are frequently asymptomatic to the infection, though they can present subclinical infections and suffer from recidivist bacteremia through their lives [6]. B. henselae transmission between cats occurs through arthropod vectors, mainly fleas [7] and dogs are also potential reservoirs of B. henselae [8]. The transmission route between humans and cats occurs mainly through the dermic inoculation of the bacteria through scratches or bites of contaminated cats [9]. However, other studies show B. henselae infections of immunocompetent humans through the possible inoculation by ticks and spider bites [10,11]. In the bloodstream, the bacteria invade erythrocytes, where they persist intracellularly, causing erythrocytic and vascular endothelial alterations [12].

B. henselae infection in humans causes multiple clinical signs, such as dermatic, cardiovascular, lymphatic, hepatic and nervous system diseases [13]. CSD is characterized by circumscribed regional lymphadenopathy in the inoculation site and previous studies revealed an immune-dependent pattern of clinical manifestations of CSD, mainly in immunosuppressed humans [14].

Although the bacteria present a high veterinary, medical, and zoonotic importance, there are currently few studies about the genomic profile of the species and the few existing ones aim at characterizing specific genes and analyzing the recombination and mutation rates [15-18]. Previous studies revealed the possibility of a horizontal gene transfer of adaptability and virulence genes in the species from the genus Bartonella [19]. Altogether, the fact that B. henselae is distributed worldwide and considered as a neglected zoonotic pathogen, opens doors for new genome plasticity studies of the species. In this work, we performed phylogenomics, pangenomics, and genome plasticity analyses to find the possible epidemiological relationships between the strains, the conserved and variable subsets of genes and the pathogenicity and antibiotic resistance islands that may be involved in the pathogenesis process.

Material and Methods

Genomes

We retrieved 24 genomes of different strains of B. henselae publicly available at the GenBank on NCBI (https://www.ncbi.nlm. nih.gov). The genomes were retrieved through FTP in .gbff, .fna and .faa file formats. The genome of Bartonella apis was used as a nonpathogenic reference where applicable.

Phylogenomic analyses

The software Gegenees was used to compare the percentage of similarity between the 24 genomes [20]. The .fna files containing the complete and draft genomes of B. henselae were imported to Gegenees and analyzed using the default parameters. Briefly, Gegenees fragments all genomes in pre-defined sizes; performs similarity analyses using BLASTn among the genomes to identify the commonly shared regions; and, finally, the software creates a heatmap with the percentage of similarity among the genomes. Here, we used a fragmentation size of 200 bp and a similarity threshold of 40%.

The similarity matrix generated by Gegenees was then exported in .nexus format and further used in the software SplitsTree [21] to create a neighbour-joining phylogenetic tree based on the similarity among the strains using the option Equal Angle.

For a greater resolution on the phylogenomic tree analysis, a whole genome Multilocus sequence typing (wgMLST) analysis [22] was performed via the online software PGAdb-builder [23]. The 24 genomes of B. henselae, in the .fasta format, were compared with the PGAdb profile through BLASTn, using the module Build_PGAdb. Next, the PGAdb profile of the genomes, in the .scheme format, was used to build a wgMLST tree using the module Build_wgMLSTtree, with a 90% coverage and 90% identity filter. After the BLASTn analyses, the output file .newick was exported and used as input file in the software MEGA7 [24] for the reconstruction of the phylogenomic tree.

Gene Synteny

The gene synteny analysis was performed using the software Mauve [25]. Briefly, Mauve fragments the genomes in pre-defined sizes, creates Locally Collinear Blocks (LCBs) using the sequence alignments and exports the results as a figure where the rearrangement events are represented. Here, we imported the .fna files from the 24 strains and used the Houston-1 strain as reference along with the “progressiveMauve” algorithm.

Pan-genomic analyses

The prediction of orthologous genes was performed with the software OrthoFinder [26] and further classified with the use of an in-house script in: core genome, containing only genes that are commonly shared by all strains; shared genome, containing genes that are shared by two or more strains, but not all; and singletons, with strain-specific genes. Briefly, OrthoFinder used the .faa files with the amino acid sequences of all Coding DNA Sequences (CDSs) in each genome to perform an all-vs-all BLASTp analysis. The sequences were then grouped using the Markov Clustering (MCL) algorithm to determine the orthologous genes [27].

Pan-genome and subsets development

The pan-genome, core genome and singletons development were calculated based on the mean values of the permutations of all genomes using the method described in Soares et al. [28]. The final curves were then fitted using an in-house script to estimate the fixed parameters for Heap’s Law (pan-genome analyses) and leastsquares fit of the exponential regression decay (core-genome and singletons). The extrapolations of the pan-genomes from the different datasets were calculated by curve fitting based on Heap’s Law with the formula n = κ *N^-a, where n is the expected number of genes for a given number of genomes, N is the number of genomes, and the other terms are constants defined to fit the curve. The extrapolations of the core genomes and singletons for all datasets were calculated by curve fitting based on least-squares fit of the exponential regression decay with the formula n = κ*exp[-x/τ]+τg(0), where n is the expected subset of genes for a given number of genomes, x is the number of genomes, exp is Euler’s number, and the other terms are constants defined to fit the curve.

We used 7 subsets of genes for all the analyses: the complete genome dataset; 2 subsets with genomes isolated from different hosts (cats and humans); two subsets of genomes from different locations (France and USA); and, two subsets containing the genomes of the two groups identified on Gegenees and PGAdb (genomes 2-8 and 9-25 on Gegenees). The group 2-8 includes the strains: A242, A244, A121, A112, A233, U4 and BM1374165, while the group 9-25 consists of the strains: A235, JK41, JK42, Zeus, JK53, A76, A74, BM1374163, MVT02, A20, A71, F1, FDAARGOS175, Houston-I, Houston-1, JK51, and JK50.

Classification of CDSs of the pan-genome subsets into the cluster of orthologous groups

The subsets of the pan-genome (core, shared and singletons) were classified according to the functional categories of the cluster of orthologous groups (COG) into 1 – Information storage and processing; 2 – signaling and cellular processes; 3 – metabolism; and, 4 – poorly characterized. To perform this classification, the CDSs of the subsets were BLAST aligned against the myva database of COG, using an e-value of 1e-6, and the result was crosschecked with the WHOG information from COG [29].

Horizontally acquired regions

The software GIPSy (Genomic Island Prediction Software) [30,31] was used to perform the prediction of Pathogenicity Islands (PAIs) and Resistance Islands (RIs) in the genomes of B. henselae. Briefly, GIPSy is a multi-pronged tool that identifies the most common features of the pathogenicity islands, such as genomic signature deviation, tRNA flanking genes, transposases and a high concentration of virulence factors. We used .gbk files and default parameters in our analyses. The genome of Bartonella apis [32] was used as the non-pathogenic reference.

The phage regions, i.e. regions harboring phage sequences, were predicted and annotated by the software PHASTER [33] using the .fasta file of the Houston-1 strain. PHASTER performs gene prediction in bacterial genome followed by BLAST analyses against a customized phage database on NCBI and prophage database developed by Srividhya et al. [34,35], which contains phage genes. PAI and phage regions were plotted in comparative circular genome representations on the software BRIG (BLAST Ring Image Generator) [36].

We used the software BRIG to create a circular genome mapping using different genomes as references. The genomic sequences in. fasta format were used on BRIG, which performed the comparative analysis using BLAST to compare the strains of B. henselae with the reference Houston-1 strain. Furthermore, the species B. apis was also added to the analysis. The results were plotted on a circular comparative map. The image generated is composed of rings, where each genome is represented by a ring of a given color. The intensity of such color is related to the degree of similarity with the reference genome. Deleted regions of the genomes, compared with the reference, are represented by empty regions (blanks). The coordinates of the islands predicted by GIPSy and the phages predicted by PHASTER were added to BRIG’s circular map for visualization of the genomic plasticity events. Two additional plasticity analyses were performed using, as references, strains isolated from cats and humans (BM1374165 and Houston-I, respectively). Both were selected as references because they are strains with the largest complete genomes among the 24 available genomes for each host group.

The coordinates of the pathogenicity islands predicted by GIPSy were analyzed at the ACT (Artemis Comparison Tool) software [37] and exported to a list of the genes present in each region. The tool PATRIC [38] was used to find information on the virulence genes found in the analysis. The tool uses BLASTp to compare the genes found on the islands against a database that integrates the Virulence Factor Database (VFDB) [39], MvirDB [40] and manually annotated virulence genes, mainly of the genus Mycobacterium, Shigella, Salmonella, Escherichia, Listeria and Bartonella [38].

Results

General Features

Abstract

Abstract

Abstract

References

1. Who. Antiretroviral Drugs For Treating Pregnant Women And Preventing Hiv Infection In Infants Recommendations for a public health approach 2010 version. 2010.
2. Aids G. Global Aids Up. 2016.
3. Tabod T. A logical step or an alienating move ? Mandatory HIV testing in Cameroon : 2016.
4. Johansson Ka, Pedersen Kb, Andersson A-K. Hiv Testing of Pregnant Women: an Ethical Analysis. Dev World Bioeth [Internet]. 2011; 11: 109-19.
5. Armstrong R. Mandatory HIV testing in pregnancy: Is there ever a time? Dev World Bioeth. 2008; 8: 1-10.
6. Kongnyuy EJ. HIV testing during pregnancy in sub-Saharan Africa. Expert Rev Obstet Gynecol. 2009; 4: 359-62.
7. Bain LE, Dierickx K, Hens K. Ethical issues surrounding the provider initiated opt - Out prenatal HIV screening practice in Sub - Saharan Africa: A literature review. BMC Med Ethics [Internet]. BMC Medical Ethics. 2015; 16: 1-12.
8. Bain LE, Tarkang EE. Intrapartum HIV screening revisited. Which approach is appropriate for Sub-Saharan Africa today: opt-in, opt-out or both? J AIDS Clin Res [Internet]. 2014; 5: 363.
9. Tudor Car L, Brusamento S, Elmoniry H, van Velthoven MHMMT, Pape UJ, Welch V, et al. The uptake of integrated perinatal prevention of mother-tochild HIV transmission programs in low- and middle-income countries: a systematic review. PLoS One [Internet]. 2013; 8: e56550.
10. Shetty AK, Marangwanda C, Stranix-Chibanda L, Chandisarewa W, Chirapa E, Mahomva A, et al. The feasibility of preventing mother-to-child transmission of HIV using peer counselors in Zimbabwe. AIDS Res Ther. 2008; 5: 1-8.
11. Kourtis AP, Lee FK, Abrams EJ, Jamieson DJ, Bulterys M. Mother-to-child transmission of HIV-1: timing and implications for prevention. Lancet Infect Dis. 2006; 6: 726-32.
12. WHO. Guidelines on HIV self-testing and partner notification: supplement to consolidated guidelines on HIV testing services. 2016; 1-104.
13. Lester P, Patridge JC, Chesney MA CM. The consequences of a positive prenatal HIV antibody test for women. J Acquir Immune Defic Syndr Hum Retrovirol. 1995; 10: 341-9.
14. Temmerman M, Ndinya-Achola J, Ambani J PP. The right not to know HIVtest results. Lancet. 1995; 345: 969-70.
15. Weiser SD, Heisler M, Leiter K, Percy-De Korte F, Tlou S, DeMonner S, et al. Routine HIV testing in Botswana: A population-based study on attitudes, practices, and human rights concerns. PLoS Med. 2006; 3: 1013-22.
16. Gruskin S, Ahmed S FLP. Provider-intiated HIV testing and counseling in health facilities – what does this mean for the health and human rights of pregnant women? Dev World Bioeth. 2008; 8: 23-32.
17. WHO. Guidance in Provider-initiated HIV Testing and Counselling in Health facilities. Geneva. 2007.
18. Siemieniuk RAC, Jadavji T GM. Limitations of opt-out HIV screening and mother–child HIV transmission. Am J Public Health. 2010; 100: 388-9.
19. Siegfried N, van der Merwe L, Brocklehurst P ST. Antiretrovirals for reducing the risk of mother-to- child transmission of HIV infection. Cochrane Database Syst Rev. Cochrane Database. 2011.
20. H. M. Prevention of Mother-To-Child Transmission (PMTCT) Briefing Note. Clin Obs Gynecol. 2001; 44: 210-25.
21. Groves A, Wyckoff S. HIV Testing During Pregnancy A Literature and Policy Review. Violence Against Women. 2008; 20080916.
22. AM B. The ethics of mandatory versus voluntary HIV testing of pregnant women. Curr Surg. 2000; 57: 165-8.
23. Armstrong R. Mandatory HIV testing in pregnancy: Is there ever a time? Dev World Bioeth. 2008; 8: 1471-8847.
24. Bayer R, Fairchild AL. Changing the paradigm for HIV testing—the end of exceptionalism. N Engl J Med. 2006; 355: 647-9.
25. Borry P, Schotsmans P, Dierickx K. What is the role of empirical research in bioethical reflection and decision- making ? An ethical analysis. Med Heal Care Philos. 2004; 7: 41-53.
26. Flamm AL, Kodish E. Eempirical bioethics research is a winner, but bioethics mission creep is a false alarm. J Clin Ethics. 2014; 25: 189-93.
27. Gostlin L, Mann JM. Towards the Development of a Human Rights Impact Assessment for the Formulation and Evaluation of Public Health Policies. Heal {&} Hum Rights J. 1994; 1: 58-80.
28. Gruskin S. Bangkok 2004. Current issues and concerns in HIV testing: a health and human rights approach. HIV AIDS Policy Law Rev. 2004; 9: 99- 103.