Metabolic phenotype of clinical and environmental Mycobacterium avium subsp. hominissuis isolates

Bacterial isolates and BIOLOG phenotype microarray

We analyzed five clinical and five environmental MAH isolates (Table 1).

We performed the BIOLOG Phenotype Microarray™ (BIOLOG, Hayward, CA, USA) according to the manufacturer’s recommendations. The technology is based on the measurement of bacterial respiration, which produces NADH (Bochner, Gadzinski & Panomitros, 2001). If bacteria are able to metabolize a specific substrate, electrons from NADH reduce a tetrazolium dye in an irreversible reaction generating a purple color in the PM plate wells. This color change is measured and recorded every 15 min by the reporter instrument OmniLog™ (BIOLOG, Hayward, CA, USA), generating a kinetic response curve for each well (Bochner, 2003; Bochner, 2009).

The ten MAH isolates were tested with the 96-wells plates PM1 to PM4, containing 190 carbon (PM1 and PM2), 95 nitrogen (PM3), 59 phosphorous (PM4) and 35 sulfur (PM4) substrates. The PM plates 1, 2 and 3 include one negative control well, in which bacteria are tested without any substrate. The PM4 plate includes two negative control wells, one for the phosphorus and one for the sulfur substrates. All isolates were tested three times. Briefly, we cultivated each MAH isolate in 30 ml of 7H10 Middlebrook medium supplemented with 10% modified ADC-enrichment (2% of glucose, 5% of BSA, 0,85% of NaCl) until an OD_{600 nm} of 0.3–0.6 was achieved (mid-logarithmic phase of growth). The use of liquid cultures in place of agar reduces bacterial clumping. Bacterial cultures were harvested by centrifugation for 10 min at 4,000 g and pellets were re-suspended in 10 ml of distilled water. Bacterial cells were starved for one night in water at room temperature to minimize false positive reactions due to nutrient accumulation in MAH cells and to ensure the use of the substrates provided by the PM plates. The following day the cells were centrifuged and re-suspended using a sterile stick in tubes containing 10 ml of GN/GP-IF-0a (BIOLOG inoculating fluid), 120 µl of 100× BIOLOG Redox Dye Mix G and 1 ml of the appropriate additive (Table 2), until 85% transmittance was reached as measured using the turbidimeter provided by BIOLOG. In order to reduce bacterial clumping, the sterile stick used for inoculation was ground against the wall of the tube. A volume of 100 μl of this final suspension was added to each of the 96 wells of the PM plates. The PM plates were then sealed to avoid drying and incubated at 37 °C in the OmniLog^® (BIOLOG, Hayward, CA, USA) incubator reader for 8 days.

As recommended by BIOLOG, we tested plates PM1 to PM4 using the same assay protocol but without the addition of bacteria in order to identify wells with abiotic dye reduction, which can generate false positive results.

Analysis of BIOLOG phenotype microarray data

The raw kinetic data were exported as CSV files using OmniLog PM file Management/kinetic Analysis module (Bochner, 2003; Khatri et al., 2013). Differences in the metabolization of the different substrates by the ten MAH isolates were investigated by analyzing the maximum height of the bacterial respiration curves (parameter A) using the R-package opm (Vaas et al., 2013). To allow comparisons across plates processed in different experimental runs, the A parameters were normalized by subtracting the well mean of the negative control (Vaas et al., 2013). Furthermore, the A parameters of the triplicates were combined by calculating the mean and discretized into “positive,” “moderate” and “negative” metabolization using the method “discrete” within the R-opm package. Substrates differentiating the isolates from each other were visualized as a heatmap generated using the R-packages heatmap.plus with the Euclidean algorithm. The heatmap displays the utilization of each substrate with a color key: yellow for strong positive metabolization, green for moderate metabolization and blue for no metabolization.

Analysis of metabolic pathways

The metabolic pathways of the two substrates of interest butyric and propionic acid have been further analyzed. Specifically, we extracted all sequences of the genes known to be associated with the pathways related to butyric and propionic acid from the KEGG pathway database (Kanehisa et al., 2016). We extracted the genes from all the M. avium subspecies (n = 8) present in the KEGG pathway database, namely: M. avium subsp. paratuberculosis K-10, M. avium subsp. paratuberculosis MAP4, M. avium subsp. paratuberculosis E1, M. avium subsp. paratuberculosis E93, M. avium subsp. avium DJO-44271, M. avium subsp. avium 2285 (R), M. avium subsp. avium 2285 (S) and the M. avium 104. The redundant genes have been excluded. Then we screened all such genes in genomes of our ten MAH isolates by performing a Custom BLAST analysis using Geneious version 9 (Kearse et al., 2012). The parameters for the screening that we used to determine if a gene was present or not were: sequence identity ≥90%, sequence coverage ≥90%, e value ≤0.01.

In addition, we analyzed the number of Single Nucleotide polymorphisms (SNP)s (both synonymous and nonsynonymous) in the sequence of the genes detected in our MAH isolates. For each gene, we also constructed a phylogenetic tree using the nucleotide sequences to determine whether any SNP was associated with clinical or environmental source of the isolates based on the Tamura–Nei model using Geneious version 9.

Statistical analyses

We generated two groups, one with data from all clinical isolates and the other with data from all environmental isolates. Statistical differences between clinical and environmental isolates in the metabolization of butyric acid and propionic acid were evaluated by means of 95% family-wise comparison of group means (Tukey contrast test) of the parameter A on specific wells using the function “opm_mcp” within the opm R-package. A p value less than 0.05 was considered to be statistically significant.

Whole genome sequencing of MAH isolates

Genomic DNAs were extracted from the MAH isolates as described previously (Lewin et al., 2003). Whole genome sequencing (WGS) was performed using Illumina MiSeq 300 bp paired-end sequencing, yielding a coverage that exceeded 100×. The NGS QC tool kit was used to assess the quality of the data reads, which was set as reads with a minimum of 70% of bases having a phred score greater than 20 (Patel & Jain, 2012). De novo assembly of the resulting reads into multiple contigs was performed using CLC Genomics Workbench 8.0 (CLC bio, Aarhus, Denmark) and contigs annotation was done using RAST (Aziz et al., 2008).

Determination of the maximum common genome and of the accessory genome

We determined the maximum common genome (MCG), comprising those genes present in all of the ten MAH genomes, as reported previously (Von Mentzer et al., 2014). All these genes were then extracted from all genomes, concatenated and aligned. The resulting alignment was used to generate a clustering tree using RAxML 8.1 (Stamatakis, 2014).

For determination of the accessory genome we applied the PanGenome Pipeline –Roary. After determination of the accessory genome of the ten MAH genomes and its distribution within them, we separated those genes that are exclusively present only in either the environmental strains or the clinical strains (Page et al., 2015).

Substrate utilization of the ten MAH isolates

We tested the capability of our ten MAH isolates to metabolize 379 different substrates. In total, 334/379 (88.1%) substrates were negative for all of the isolates (see Table S1). A total of 23/379 (6.1%) substrates caused abiotic reactions and were excluded from further analysis. A list of false-positive substrates is shown in the Table S2. The kinetic curves corresponding to the control plates PM1 to PM4 tested without bacteria are presented in the Fig. S1.

Only two carbon substrates, the fatty acid derivatives Tween 20 and Tween 40 were strongly positive for all of the ten MAH isolates. The kinetic curves for these substrates reached 250 Omnilog units, amongst the highest values recorded in our analysis (see Fig. S2 for all kinetic curves of the ten MAH isolates). The opm analysis revealed that a total of 20/379 (5.3%) substrates were metabolized differently among the MAH isolates (Table 3). We therefore carried out further analysis using only these substrates. The majority of these 20 substrates were carbon substrates, 15/20 (75.0%), followed by 3 nitrogen and 2 phosphorous substrates. The heatmap in Fig. 1 shows the utilization of these 20 substrates among the ten MAH isolates. The isolates are grouped according to their substrate utilization. Isolates utilizing similar substrates appear to cluster together.

Two major clusters, each composed of five isolates, could be observed. One was rich in environmental isolates (4/5) and the other was rich in clinical isolates (4/5). The substrates predominantly contributing to this clustering were butyric acid and propionic acid and indeed, the Tukey’s test revealed that environmental isolates metabolized more strongly butyric acid (p = 0.0209) and propionic acid (p = 0.00307) than clinical isolates with statistical significance.

Metabolic pathways analysis

The propionic and butyric acid are involved in three and one pathway, respectively (Table 3). A total of 151 genes have been identified in the KEGG database associated with all these pathways (Kanehisa et al., 2016). In Table S3 we reported the distribution and SNPs analysis of those genes in the MAH genomes. Of the 151 genes, 134 (88.7%) are present in all the ten MAH. The median gene length was 1,099 bp (range 318–2,253), whereas the median number of SNPs per gene is 17 (range 1–147). The phylogenetic analysis revealed that none of the SNPs could be associated with the group of the clinical or the group of environmental MAH isolates (see Table S4). In the propanoate pathway four operons have been identified: fadAB associated with the β-oxidation of several fatty acids (DiRusso, 1990), ech8-9 encoding for hydrogenases that play a role in energy conversion (Sant’Anna et al., 2015), sucCD responsible for the succinate metabolism (Cerdeno-Tarraga et al., 2003) and mutAB involved in the methylmalonate pathway (Schoenwolf & Alvarez, 1989). In the nicotinate pathway there are two operons: the pntAA-AB-B responsible for the transhydrogenation between NADH and NADP (Anderlund et al., 1999) and nadABC involved in the biosynthesis of NAD+ (Vilcheze et al., 2010). In the degradation of aromatic compounds pathway we identified the pcaHGB operon involved in the β-ketoadipate pathway (Harwood & Parales, 1996). In the butanoate pathway we identified the fadAB and ech8-9operons, the sdhCDAB encoding for the succinate dehydrogenase complex involved in the fatty acid metabolism (Nam et al., 2005) and the ilvBN responsible for the acetolactate synthesis, a precursor of several amino acids (Keilhauer, Eggeling & Sahm, 1993).

Clustering analysis and determination of the accessory genome

The WGS of the two reference strains MAH 104 and MAH TH135 were already in the GenBank database and we submitted the remaining genomes at DDBJ/EMBL/GenBank under the BioProject Number PRJNA299461. The MCG, the maximum number of genes shared by all ten MAH isolates was 1,658, the alignment of which spanned 1.378 Mbp. The clustering analysis of the ten MAH isolates is shown in Fig. 2. By comparing the genetic clustering obtained by WGS with the phenotypic clustering obtained through BIOLOG PM we observed slight differences. For examples, the isolates MAH E-96-2 and MAH E-82-7, which share identical metabolic profiles, were genetically more distant from each other. Interestingly, at the genetic level there was no obvious clustering between the group of clinical and the group of environmental isolates.

The accessory genome is constituted by 4,067 genes. A total of 1,688 genes were specific for the group of clinical isolates (Table S5). On the other hand, 698 genes were specific for the group of environmental isolates (Table S6). We found no genes that were present in all the clinical and absent in all the environmental isolates, and vice-versa. Among the most abundant specific genes of the two groups of isolates, there were no known genes associated with the pathways which involved butyric and propionic acid. However, genes annotated as hypothetical proteins represented the most abundant specific genes of the two groups