Genetic characterization of Bacillus anthracis strains circulating in Italy from 1972 to 2018

Summary

In Italy anthrax is an endemic disease, with a few outbreaks occurring almost every year. We surveyed 234 B. anthracis strains from animals (n = 196), humans (n = 3) and the environment (n = 35) isolated during Italian outbreaks in the years 1972–2018. Despite the considerable genetic homogeneity of B. anthracis, the strains were effectively differentiated using canonical single nucleotide polymorphisms (CanSNPs) assay and multiple-locus variable-number tandem repeat analysis (MLVA). The phylogenetic identity was determined through the characterization of 14 CanSNPs. In addition, a subsequent 31-loci MLVA assay was also used to further discriminate B. anthracis genotypes into subgroups. The analysis of 14 CanSNPs allowed for the identification of four main lineages: A.Br.011/009, A.Br.008/011 (respectively belonging to A.Br.008/009 sublineage, also known Trans-Eurasian or TEA group), A.Br.005/006 and B.Br.CNEVA. A.Br.011/009, the most common subgroup of lineage A, is the major genotype of B. anthracis in Italy. The MLVA analysis revealed the presence of 55 different genotypes in Italy. Most of the genotypes are genetically very similar, supporting the hypothesis that all strains evolved from a local common ancestral strain, except for two genotypes representing the branch A.Br.005/006 and B.Br.CNEVA. The genotyping analysis applied in this study remains a very valuable tool for studying the diversity, evolution, and molecular epidemiology of B. anthracis.

Keywords:

Polymerase chain reaction – Phylogeography – Veterinary diseases – Italian people – Italy – Bacterial spores – Anthrax – Bacillus anthracis

Introduction

Anthrax is a non-contagious zoonotic disease affecting a broad range of animal species including humans. Bacillus anthracis, the etiological agent of anthrax, forms highly resistant spores that can to persist in the environment for several decades [1]. Domestic and wild ruminants are species most susceptible to anthrax [2]. Animals are infected during grazing in areas contaminated with anthrax spores, while humans can contract the disease by contact with anthrax-infected animals or anthrax-contaminated animal products. Most frequently this involves employment in specific high risk occupation; such a farmer, butcher, tanner, wool carder, shearer and veterinarian. Exposure most commonly occurs during the skinning and butchering of cattle that are sick or dead because of anthrax [3]. Three forms of anthrax occur in humans, depending on exposure type: cutaneous (which is usually non-fatal), gastrointestinal, and inhalational, both of which can be potentially fatal [4]. Recently, a fourth form of the disease was reported in drug users who inject drugs contaminated with anthrax spores [5]. Further, since it is relatively easy and inexpensive to obtain B. anthracis, the bacterium is one of the preferred pathogenic agents for use as a bacteriological weapon in bio-terrorist attacks [6]. In Italy, anthrax is typically a sporadic disease, particularly occurring during the summer (with a few exceptions) in the central and southern regions, and in the major islands, where it almost exclusively affects animals at pasture [7]. Between 1972 and 2018, approximately 200 outbreaks of animal anthrax were recorded (unpublished data). Very rarely, anthrax infection takes the form of an epidemic-like disease, characterized by outbreaks in limited areas involving a great number of animals. In Italy, two major epidemic-like anthrax outbreaks have been reported: one during the summer of 2004 in Basilicata, and one during the summer of 2011, in an area between Basilicata and Campania [8, 9]. Molecular tools, such as the canonical SNPs assay (CanSNPs) and multiple-locus variable-number tandem repeat analysis (MLVA) are highly effective for differentiating B. anthracis strains. The overall goal of this study was to utilize SNP analysis to establish the phylogenetic relationship between the B. anthracis strains examined, and further discriminate them in the context of the MLVA assay, in order to examine correlations among the B. anthracis isolates associated with the Italian anthrax outbreaks and to assess genetic diversity at regional and broader scales.

Materials and methods

Ethics statement

The animal and environmental strains used in the current study were isolated at the Anthrax Reference Institute of Italy (Ce.R.N.A.), a public laboratory mandated by the Italian Ministry of Health to confirm diagnosis of all animal anthrax cases in Italy. During outbreaks, samples were collected by the veterinary services of the Ministry of Health. Specific permission for soil sampling was not required. B. anthracis DNAs from anthrax human cutaneous cases were also included in the current study; two came from the “San Carlo” Hospital, Department of Infectious Disease, Potenza, Italy, and one from the “L. Spallanzani” National Institute for Infectious Disease, Rome, Italy [10].

Bacterial strains

A collection of 234 B. anthracis strains, including 196 strains isolated from animal and 35 from the environment, isolated during Italian anthrax outbreaks in the years 1972–2018, were analyzed in the current study (Table 1). Furthermore, 3 B. anthracis DNAs from anthrax human cutaneous cases were also analyzed.

DNA extraction

B. anthracis strains were seeded on 5% sheep blood agar plates and then incubated at +37°C for 24 h. Bacterial DNA was extracted using the DNAeasy Blood and Tissue kit (Qiagen, Hilden, Germany), following the protocol for Gram-positive bacteria. All manipulations of B. anthracis strains were performed in a biosafety level 3 laboratory at the Experimental Zooprophylactic Institute of Apulia and Basilicata Regions in a class II type A 2 biosafety cabinet.

Real-time polymerase chain reaction (PCR) assay

Molecular identification of B. anthracis was performed using qualitative real-time PCR. The method is based on the amplification of specific DNA sequences using three pairs of specific primers [11] as follows: R1/R2 primers, specific for the BA813 gene located on the B. anthracis chromosome; PAG 23/24 primers, specific for the protective antigen (PA) gene located on the virulence plasmid pXO1; and CAP 57/58 primers, specific for the capsule (CAP) gene located on the virulence plasmid pXO2. Each 20 μl reaction mixture contained 1x Sso Advanced TM SYBR® Green Supermix (BIORAD), 300 nM each forward and reverse primer, and approximately 10 ng DNA template. The amplification was performed using the CFX Connect Real Time PCR Detection System (BIORAD). A melting curve was generated at 0.5°C increments between 65°C and 95°C, and was analyzed by CFX Manager TM Software, Version 3.0 (BIORAD).

CanSNP analysis

CanSNP profiles were obtained using 13 allelic discrimination assays involving specific oligonucleotides and probes, as described by Van Ert et al. [12]. Each 10 μl reaction mixture contained 1x TaqMan Genotyping Master Mix (Applied Biosystems, Foster City, CA, USA), 250 nM probe, 600 nM each of forward and reverse primer, and approximately 10 ng DNA template. For all assays, the thermal cycling parameters used were as follows: 10 min at 95°C, followed by 40 cycles of 15 s at 95°C and 1 min at 60°C. Endpoint fluorescent data were acquired by using the ABI 7900HT instrument. The CanSNPs data were compared with the data for 12 recognized sublineage or subgroups. The 14th SNP was detected using a High Resolution Melting (HRM) assay for a specific A.Br.011 CanSNP [13,14]. Position 2,552,486, based on the Ames Ancestor genome (NC_007530.2), was analyzed. Amplification was performed using the CFX Connect Real Time PCR Detection System (BIORAD) and Precision Melt Supermix (BIORAD). The reaction mixture contained 0.2 μM of each primer and 1x Precision Melt Supermix (BIORAD) in a 20 μl final volume. The following cycling parameters were used: 2 min at 95°C, were followed by 35 cycles of 10 s at 95°C and 30 s at 60°C. The samples were then heated to 95°C for 30 s, cooled down to 60°C over 1 min, and then heated from 65°C to 95°C at a rate of 0.5°C/s. High-resolution melting data were analyzed using Precision Melt Analysis Software (BIORAD).

31-loci MLVA analysis

For the 31-marker MLVA, 5' fluorescently labeled oligonucleotides (6-FAM, VIC, NED and PET), specifically selected for variable number tandem repeats (VNTR) analysis were used. Twenty-seven chromosomal VNTR loci (vrrA, vrrB1, cg3, vrrB2, vntr19, vrrC1, vrrC2, vntr32, vntr12, vntr35, vntr23, bams03, bams05, bams13, bams15, bams21, bams22, bams23, bams24, bams25, bams28, bams30, bams31, bams34, bams44, bams51 and bams53) and four plasmid loci (vntr16, vntr17, pxO1 and pxO2) [12, 15–18] were analyzed. The MLVA assay involved preparation of two singleplex and nine multiplex reactions, in a final volume of 15 μl. Each reaction mixture contained the following: 1x PCR reaction buffer (Qiagen, Hilden, Germany); 3 mM MgCl₂, 0.2 mM for each dNTP; 1 U Hot Star Plus Taq DNA polymerase (Qiagen, Hilden, Germany), the appropriate concentration of each primer (as described in Table 2); and approximately 10 ng DNA template.

The following PCR cycling program was used for the two singleplex reactions and for multiplex reactions 1 and 2: 5 min at 95°C; followed by 35 cycles of 30 s at 94°C, 30 s at 60°C, and 30 s at 72°C, with a final step of 5 min at 72°C. The following amplification program was used for multiplex reactions 3: 5 min at 95°C, followed by 35 cycles of 30 s at 94°C, 30 s at 54°C, 45 s at 72°C, and 5 min at 72°C. The following amplification program was used for multiplex reaction 4: 5 min at 95°C, followed by 35 cycles of 30 s at 94°C, 45 s at 56°C, 1 min at 72°C, and 5 min at 72°. The following amplification program was used for multiplex reaction 5: 5 min at 95°C, followed by 35 cycles of 30 s at 94°C, 45 s at 59°, 1 min at 72°C, and 5 min at 72°C. The following amplification program was used for multiplex reactions 6 to 9: 5 min at 94°C, followed by 35 cycles of 1 min at 94°C, 90 s at 60°, 90 s at 72°C, and 15 min at 72°C.

Automated genotype analysis

The MLVA PCR products were diluted 1:80 and analyzed by capillary electrophoresis using the ABI Prism 3130 genetic analyzer (Applied Biosystems) and 0.25 μl GeneScan 1200, and were sized by using Gene Mapper 4.0 (Applied Biosystems Inc.). Assignment of the sizes and corresponding repeating unit numbers for each locus was performed using the following strains as references: Ames Ancestor (NCBI Reference Sequence: NC_007530.2, chromosome), pXO1 (NCBI Reference Sequence: NC_007322.2, plasmid), and pXO2 (NCBI Reference Sequence: NC_007323.2, plasmid). Data were analyzed using conventional values proposed in the updated version of the 2016 Bacillus anthracis MLVA database, available at MLVAbank (http://mlva.u-psud.fr/). A phylogram was derived by clustering with the unweighted pair group method with arithmetic means (UPGMA), using ‘categorical’ character table values. All markers were given equal weight, irrespective of the number of repeats.

The discriminatory ability of the MLVA technique was determined by calculating the discriminatory index (D) for 234 typed strains. The discriminatory power of each VNTR was estimated by the number of alleles detected and the allele diversity using BioNumerics 7.6 software (Applied Maths, Belgium) [19].

Results

Real Time PCR, CanSNPs and MLVA analysis of anthrax strains

All the analyzed strains tested positive after the PCR amplification of chromosomal, plasmid pXO1 (toxins coding) and pXO2 (capsule formation) targets. The analysis of 13 classical CanSNPs revealed that 231 analyzed strains belonged to the sublineage A.Br. 008/009, also known as Trans-Eurasian (TEA) group. The TEA group was established in southern and eastern Europe and represents the dominant subgroup in Italy, Bulgaria, Hungary and Albania [7, 12, 20–22]. The analysis of an additional 14th CanSNP (A.Br.011), recently allowed for the differentiation of the A.Br.008/009 group into 2 subgroups. Accordingly, the analysis of the 14th CanSNP in the current study revealed that 207 of the 231 strains belonged to the main sub-lineage A.Br.011/009, while 24 strains belonged to the sublineage A.Br.008/011. All but one strain belonging to the latter sublineage were isolated in Sicily; one strain was isolated in Umbria. Further, one strain isolated in Veneto belonged to the main lineage A, sublineage A.Br.005/006, while two other strains, one from Veneto and one from Trentino, belonged to the main lineage B, sublineage B.Br.CNEVA.

MLVA based on the analysis of 31 VNTRs revealed 55 different genotypes, as shown in S1 Table, distributed in the Italian regions named GT-1 to GT-55, accordingly (Fig 1). The GT-14 genotype was the most common and was represented by 34 strains, mostly from Basilicata, Apulia, and Calabria. The GT-31 genotype was represented by 19 isolates: 16 from Tuscany, two from Apulia and one from Sardinia. The GT-26 and GT-27 genotypes were only isolated in the Basilicata and Campania regions. Other genotypes were characteristic for single regions, as showed in Table 3.

The number of different alleles ranged from 1 for bams21 and bams25 to 10 for bams15. Highest allelic diversities measured by Shannon Diversity Index (0.40632) was observed for the locus bams15 (Table 4). The relationship among the strains based on MLVA results is represented in Fig 2.

Discussion

Bacillus anthracis is clonal in nature and often exhibits a high degree of genetic homogeneity due to the fact that is has a single stranded chromosome and reproduces asexually. This characteristic has traditionally made the discrimination of isolates for epidemiological purposes difficult. Furthermore the high survivability of spores in the soils, allowed B. anthracis to reproduce a relatively limited number of times during its evolution [23]. The 31-loci MLVA analysis carried out on 234 B. anthracis strains, isolated in Italy during the years 1972–2018, revealed the circulation of 55 B. anthracis genotypes. The performed CanSNPs analysis placed 53 of the 55 identified genotypes in a common cluster (TEA). The analysis of the classical 13 CanSNPs revealed that most of the analyzed strains (98%) belonged to the sublineage A.Br.008/009 (the TEA group), which is the most common group in Europe and Asia [15]. However, except for the genotypes of strains isolated in Umbria and some others isolated in Sicily belonging to sublineage A.Br.008/011, all strains belonged to the sublineage A.Br.011/009. Interestingly, genotype GT-54 isolated in Veneto was very different from the other characteristic Italian strains. CanSNPs analysis confirmed this observation placing this genotype in the branch A.Br.005/006. This branch is generally present in the central-southern Africa, although it was also identified in Europe [12, 24].Furthermore, genotype GT-55; B.Br.CNEVA, isolated in Veneto and Trentino is highly differentiated from most other Italian strains examined here. This genotype is widespread in Europe and found in France, Switzerland and Germany [12, 25, 26]. In Italy, the population of B. anthracis is mainly divided into two sublineages: A.Br.011/009, definitely the most common and A.Br.008/011 present only in Umbria and Sicily. Both these sublineages belong to the large TEA group (Fig 2). The TEA group A.Br.008/009 contains a B. anthracis subpopulation that is well adapted to the northern hemisphere and predominant in Europe, Russia, Kazakhstan, Caucasus and western China [12, 27]. It has also been detected in Africa [18, 28]. This group is thought to have given rise to the western north American sublineage (A.Br.WNA), which is dominant in central Canada and much of the western USA. The presence of strains belonging to sublineages A.Br.008/011 and A.Br.011/009 might represent an effect of evolution on a common ancestral strain at the territorial level. In particular, A.Br.008/011 represents a rare and deep branching sublineage, also observed in Bulgaria, France and Turkey [29]. The spread of the TEA group to Europe and Asia is postulated to be linked to animal handling along the ancient East-West commercial routes of the Silk Road [30]. In the current study, strains belonging to the B.Br.CNEVA lineage were isolated in a relatively small area of north-eastern Italy. The relatively low diversity between the two strains demonstrated in the current study is consistent with a single introduction event of the B.Br.CNEVA lineage into the country, followed by ecological establishment and progressive in situ differentiation around the Italian Alps area [21]. Consistent with this hypothesis, the Italian strains form a cluster that is distinct from the other European B.Br.CNEVA strains. Identification of one A.Br.005/006 strain in Italy could be associated with the trade exchanges dating back when city states competed for trade and commerce throughout the Mediterranean [7]. This subgroup is well represented in Africa, but rare in Europe [12]. It is therefore quite surprising that past importations of ill or dying animals or spore-infected items from Africa, the Middle East, or even Asia, did not have a greater impact on the genetic structure of B. anthracis in the region. The higher variety of B. anthracis genotypes identified in southern Italy relative to genotypes from other Italian regions may be explained by the differences in the breeding systems between northern and southern Italy. In southern Italy, many livestock farmers use extensive farming methods, which increases the chances of grazer exposure to historical spore sites and deposits. The possibility of exposure is lower in northern Italy because most livestock farmers use intensive breeding systems. Another observation from the current study was that the neighboring regions share just a few genotypes. In particular, the GT-24 genotype was present in Apulia, Basilicata and Calabria; the GT-26 and GT-27 genotypes were identified in Basilicata and Campania; and the GT-55 genotype was identified in Veneto and Trentino. Noteworthy and difficult to explain is the dislocation of genotype GT-31, identified in Apulia, Tuscany and Sardinia. These are not neighboring regions; on the contrary, they are quite far from one another. Also in this national scenario one of the explanations could be the trade of animals or animal products within the country over the years. Nevertheless, since most genotypes are exclusive to each region, it appears that Italian B. anthracis strains may be autochthonous for a single territory. Interestingly, genotypes exclusive to specific regions were detected especially in Sicily and Sardinia, probably because of low animal movements between these islands and the rest of Italy. The analysis of chromosomal and plasmid hypervariable regions using such methods as MLVA constitutes a valuable approach for studying the diversity, evolution and molecular epidemiology of B. anthracis. Therefore, MLVA is a valid method that enables the understanding of the distribution of B. anthracis within a country.

Supporting information

S1 Table [xlsx]
Allele distribution of the 55 genotypes identified using 31 VNTR analysis.