BACKGROUND
haemophilus parasuis is an important pathogen of pigs  <cit>  but is also a frequent commensal of their upper respiratory tracts.
it causes polyserositis, meningitis and septicaemia  as well as pneumonia and pleurisy  <cit> .
historically, h. parasuis caused disease in recently weaned piglets, but with the intensification of production it now affects a wide age-range  <cit> .
h. parasuis is the leading cause of mortality, alongside the prrs virus, in nursery herds in the usa, and is the third-most important bacterial pathogen affecting finisher herds  <cit> .
h. parasuis also contributes to multi-factorial porcine respiratory disease complex, the leading cause of mortality in grower-finisher pigs in the usa  <cit> .
h. parasuis is endemic in all pig farming countries and has a major economic impact on the global pig industry, with costs including those for production and stock losses, carcass disposal, and vaccination, as well as requiring antibiotic therapy .
the study of h. parasuis to date has classified the bacterium, originally using a gel immuno-diffusion assay and more commonly now an indirect haemagglutination assay, into  <dig> different serovars with non-typeable isolates also frequently found  <cit> .
the designation of serovar in h. parasuis is thought to be predominantly due to the presence of a particular polysaccharide capsule which is specific to each serovar  .
some serovars, especially  <dig>   <dig> and  <dig>  cause disease more commonly than others .
however, no absolute relationship between virulence and serovar has been found  <cit> .
the onset of disease has been linked to many factors, including physiological stress and immune status, but the current study focuses on understanding the bacterial genetic factors that might correlate with the establishment of disease.
several experiments have tried to associate the bacterium’s genotype with phenotype, but with limited success .
experimental evidence to support definitive virulence factors in h. parasuis is sparse and most studies have been limited to small numbers of isolates or serovars.
it is likely that the relationship between genotype and virulence phenotype is complex .
the first h. parasuis genome to be sequenced was in  <dig>  <cit> , with several further sequences released since then  <cit> .
these have been used to design several experiments, including searches for novel immunogenic proteins  <cit>  and virulence factors  <cit> , but these few strains may not be reflective of the wider diversity of isolates that exist in pig populations.
previous studies show that the pan-genome of a bacterial species can be vast .
recombination and lateral transfer are involved in creating the pan-genome, which can be resolved into the core and accessory genomes.
the boundaries of the core genome can be extrapolated from a small number of genomes, but it is preferable to analyse a large number of isolates for a more comprehensive definition .
an extensive amount of variation can often be found within the accessory genome of a species and it is important to investigate this source of variation to understand the adaptive potential of a given species  <cit> .
reductions in cost have allowed the whole genome sequencing of many strains of the same species, and this can provide the necessary statistical power required to find associations between genotypes and phenotypes such as virulence.
this variation is now being explored in many pathogens including campylobacter jejuni <cit> , escherichia coli <cit> , haemophilus influenzae <cit> , staphylococcus aureus <cit>  and streptococcus pneumoniae <cit> .
such studies seek associations between phenotypes and genotypes that would inform surveillance efforts, improve diagnostics, and aid understanding of the pathogen, potentially yielding new treatments of the associated diseases.
in order to disentangle the genetic and environmental bases of glässer’s disease, we undertook a pan-genome analysis of  <dig> isolates of h. parasuis with high quality clinical metadata, covering all serovars and isolates from disease- and non-disease-causing backgrounds.
we have investigated the key features of the pan-genome and assessed the relationship between the genetics of this bacterium and collated metadata, using statistical analyses of single nucleotide polymorphisms  in the core genome, and the presence of genes in the accessory genome.
we compared our analyses to the distribution of previously published virulence factors of h. parasuis within our collection of isolates.
we found that the composition of the accessory genome was a significant factor in determining whether isolates were likely to cause disease.
our findings emphasise the importance of describing the pan-genome of a bacterial species in order to understand pathogenesis, and to underpin a new generation of vaccine and diagnostic candidates for endemic bacterial pathogens.
RESULTS
epidemiology of the isolate collection
the  <dig> isolates in the collection originate from the uk , argentina , belgium , china , denmark , germany , greece , italy , japan , spain , sweden , switzerland  and the usa , with nine isolates of unknown origin.
our analysis included the reference strains of the  <dig> serovars and the published genomes of the us strain  <dig>  and the chinese strains sh <dig> and zj <dig>  <cit> .
the dates of collection ranged between 1991– <dig> for european isolates and 1955– <dig> for the reference strains  <cit> .
for the uk isolates, we collected clinical metadata including veterinary diagnostic information from necropsy, the tissue of isolation and details of any other signs of disease in the pig.
clinical metadata for the other strains included tissue of origin, clinical signs and, for a subset of strains, experimental reproduction of disease.
we used this information to categorise the isolates into clinical and non-clinical categories.
a non-clinical isolate was an isolate identified in the upper respiratory tract in the absence of any signs of respiratory or systemic disease typical of infection with h. parasuis.
a clinical isolate was an isolate identified from an animal displaying signs of h. parasuis-related disease and isolated from a respiratory or a systemic site.
we divided the clinical isolates further into respiratory and systemic, to reflect disease severity.
respiratory isolates were obtained from the lung without signs of systemic disease elsewhere in the animal.
systemic isolates were obtained from a systemic site e.g.
brain, joint or blood which are sites typically affected in glässer’s disease.
isolates that were taken from the lung of an animal that also showed signs of systemic disease were first tested in the statistical analyses as respiratory isolates and then as systemic.
this was to allow for the possibility of these isolates being from a co-infection in the same animal or where the same isolate caused disease elsewhere in the same animal.
one hundred and seventeen of the isolates had been previously serotyped  and these included examples of all fifteen known serovars.
for details of the isolate collection including disease categories and serotyping results see additional file 1: table s <dig>  analysis of the distribution of isolates by serovar and disease association  showed that the most prevalent serovar in this collection was serovar  <dig>  followed by non-typeables  and serovars  <dig> and  <dig>  which fits with global prevalence data .
our data support the previous limited studies that show that certain serovars, particularly  <dig> and  <dig>  were more strongly linked to disease  <cit> .
while serovar is still important for the formulation of vaccine strategy against this bacterium, its relationship with virulence has long been discussed but never proven  <cit> .
the fact that there are non-clinical isolates within the same serovars highlights the limitation and uncertainty of serotyping as a technique for predicting disease potential.
we expect that there are isolates amongst the non-clinical collection that will have the potential to cause clinical disease given a different host or environment, and also that some of the clinical isolates caused disease due to other environmental factors or co-infections.
the main reason for this study was to investigate the balance between the host, environment and genotypes for example serovar determining factors, predicted virulence genes and to discover new gene associations with phenotypes such as virulence.
we propose that molecular epidemiological markers that relate to severity of disease would be more useful for h. parasuis surveillance now that rapid genome sequencing of isolates is possible.
to achieve this goal we undertook pan-genome and statistical analyses with the two-fold aim of identifying genes that correlate with virulence and so could be worth investigating as markers of disease-causing potential, and with serovar to allow the design of a molecular-based serotyping technique.figure  <dig>
h. parasuis
isolate collection displayed by serovar and disease.
each serovar includes those found with cross-reactions.
nt denotes 'non-typeable’ by serotyping either due to no reaction or due to more than three reactions to the serotyping antisera.
strains without serotyping data are excluded.
the pan-genome of h. parasuis
we use the term pan-genome to refer to the set of all of the genes predicted from every isolate in the collection.
the pan-genome comprised homology groups encoding  <dig>  genes from the orthomcl analysis of the  <dig> isolates.
of these,  <dig>  genes were found within all isolates and so were designated as the core genome.
this core genome is approximately  <dig>  mb of the  <dig>  mb average genome size with  <dig>  snps identified.
the accessory genome was made up of  <dig>  homology groups that varied between isolates.
<dig> of these homology groups were classed as pseudogenes , and  <dig> % of these homology groups were identified in less than 1% of the isolates in the collection.
for details of the cog functional classification of the pan- genome see additional file  <dig> and additional file 3: figure s <dig>
synteny of the pan-genome
the syntenic pan-genome of h. parasuis is the pan-genome reordered to reflect the chromosomal location of each gene based on the sh <dig> strain  <cit> .
genes were inserted from contigs in descending size order, based on their order along each contig.
as an iterative process not all genes could be assigned to a location in the syntenic pan-genome.
therefore the syntenic pan-genome of h. parasuis represents the  <dig>  genes out of the total genes where the position could be identified in the pan-genome order, and is shown in additional file 4: figure s <dig>  several regions of variation can be seen within additional file 4: figure s <dig> with the location of many phage genes, highly variable regions and antibiotic resistance genes denoted by colouring of the heat-map.
the capsule genes did not cluster in one location in the syntenic pan-genome, which may be due to recombination in the genomes, and so the synteny of the capsule locus was established separately and can be seen in additional file 5: figure s <dig>  several phages were identified within this collection including the mu phage , a siphovirus  as well as other phage genes throughout the pan-genome.
for details of genes identified within these variable regions see additional file  <dig>  briefly, within these regions we identified many phage genes, two conjugal transfer operons, adhesins, toxins, iron-sulphur binding proteins and many proteins of unknown function.
based on the similarity of the g + c content to the average for this bacterium, many of the genes in region  <dig> appear to be from h. parasuis genes that have yet to be characterised and submitted to ncbi.
regions  <dig>   <dig> and  <dig> contained some genes of relatively high g + c content  but no hits were found when blast analysis was performed for these encoded proteins using the ncbi non-redundant protein database .
only region  <dig> contained multiple proteins with hits from nr to other bacterial species including actinobacillus, gallibacterium, klebsiella and mannhaimia.
the average g + c content of the pan-genome was  <dig> % g + c content , ranging between  <dig> - <dig> % g + c content for all isolates, which is discussed further in additional file  <dig> and additional file 6: figure s <dig>
population structure of h. parasuis
identification of two divergent clades of h. parasuis
analyses of the population structure of this isolate collection was performed using bayesian analysis of population structure   <cit>  on the core genome.
five baps populations were identified from these analyses, these represented groups within the overall population that share similar snps and indicate a common ancestry with similar amounts of admixture.
the baps plot shown in additional file 7: figure s <dig> demonstrates the within-cluster and between-cluster variation for the five baps populations based on the core genome.
while the partitioning of population structure in other bacteria has been reported before , further analysis of the causes of the partitioning usually leads to an explanation dependent on ecological separation e.g.
due to geographical separation or restriction to different hosts or tissues  <cit> .
h. parasuis has only been identified in pigs and wild boar, and the route of infection is via the nose and the upper respiratory tract and so are all able to inhabit the same niche  <cit> .
furthermore, the global nature of the pig industry means that import and export of stock occurs worldwide and so we did not expect that geographical separation would have much influence on the population structure.
we began our investigation into the population structure by checking these potential ecological barriers by comparing our clinical metadata  to the phylogenetic tree of the h. parasuis core genome built from the concatenated core genes   <cit> .
the topology of the phylogeny shows the segregation of the tree into two clades, which fits with the separation of the baps populations with baps populations  <dig> and  <dig> found in clade  <dig> and groups  <dig>   <dig> and  <dig> in clade  <dig>  we also investigated the correlation between the baps populations  and the clinical metadata .
when both the core genome and baps populations  were compared to the metadata, some separation of the serovars can be seen between the two clades, with serovars  <dig>   <dig> and  <dig> identified predominantly in clade  <dig>  serovar  <dig> is the only serovar to be found in equal proportions in both the clades.
serovar  <dig> is mostly found in clade  <dig>  for the remainder of the serovars, the number of isolates was too low to be able to restrict them to one clade.
in comparison, no influence of disease association on population structure was observed, as the clinical isolates were found throughout the tree and the baps populations .
this is contradictory to the population grouping predicted by the multi-locus sequence typing  scheme  <cit> , which showed six main subgroups on the mlst phylogeny; one predominantly clinical, one mostly non-clinical and four mixed clades.
analysis of the individual baps populations in figure  <dig>  showed that baps populations  <dig> and  <dig> contained very few non-clinical isolates .
some geographic structure can be seen in figure  <dig>  with a greater proportion of the uk isolates  in clade  <dig> which is likely to be due to the large proportion of uk isolates in the collection.
with the global nature of the pig market it is not unexpected that isolates from different geographic locations or serovars would be found in each of the baps populations.
indeed, populations  <dig> and  <dig> are very mixed; indicating only limited separation by a geographic barrier.
the baps populations are also very mixed by tissue of isolation and so no tissue tropism can be seen based on the population structure.
of the categories of metadata, only serovar shows some association with the baps populations and the differences in recombination.
therefore we conducted a second level baps analysis which split the current five baps populations into more defined populations; however the second level baps analysis still did not explain the population structure .figure  <dig>
core genome neighbor-joining tree  of
h. parasuis
.
trees are overlaid with the populations from the bayesian analysis of population structure , which represent isolates with similar rates of homologous recombination.
further metadata including disease association, serovar and country of origin are also shown.
baps populations explain the separation of the isolates on the tree into two main clades.
the baps populations defined using bayesian analysis of population structure represent isolates with similar rates of homologous recombination.
the baps populations have been compared to serovar, disease association and country of origin.
some similarity in the baps populations  <dig> and  <dig> can be seen based on serovar distribution.
very little difference can be seen between baps population when disease association is considered.
uk isolates are overrepresented in baps populations  <dig>   <dig> and  <dig>  while  <dig> and  <dig> are of mixed geographic origin.
high rates of recombination were identified in h. parasuis
analysis of the regions undergoing recombination in the core genome was performed using bayesian recombination tracker   <cit>  on the two separate clades present in figure  <dig>  we decided to perform separate analyses because the long-branch separating clades  <dig> and  <dig> would reduce the statistical power to detect recombination in either clade in a joint analysis.
the recombination analysis was also performed using gubbins  <cit> , which was performed on the alignment for all isolates.
we were able to calculate the recombination rate relative to the mutation rate  or homologous recombination rate for the individual baps populations using the output from gubbins; baps  <dig> had the highest r/m of  <dig>  followed by baps  <dig> , baps  <dig> , baps  <dig>  and finally baps  <dig> .
however, each baps population contained isolates with a range of recombination rates, with multiple isolates with a r/m of  <dig>  when the second level baps populations were considered, some clusters of isolates with r/m of  <dig> were identified but the majority of these isolates were still separated between the  <dig> groups identified in the second level analysis.
even with these isolates with low estimates of recombination rates, the r/m for each baps population is very high compared to published recombination rates in bacteria and is far higher than previously estimated for h. parasuis.
one previous study calculated the r/m for h. parasuis based on the genes involved in the mlst scheme  <cit> , and identified higher rates for the clinical isolates  compared with the non-clinical isolates .
our estimates for the rates of homologous recombination in h. parasuis far exceed the estimates based on the mlst housekeeping genes.
this may be due to purifying selection that is likely to be affecting housekeeping genes and so the estimate based on the core genome, which contains far more genetic information, should give a more informative estimate of recombination across the genome  <cit> .
the values of r/m of the two clades of h. parasuis estimated using both recombination analysis methods were within an order of magnitude.
the r/m for the first clade was  <dig> and  <dig> , while the r/m for the second clade was estimated at  <dig>  and  <dig> , using gubbins and bratnextgen respectively.
while both methods yielded comparable estimates of homologous recombination rates, the diversity in this h. parasuis population is greater than is usually included in the analysis using gubbins and so care should be taken when interpreting the results.
gubbins was designed to identify recombinations from distant sources within very clonal populations and so these r/ms may even be an under-estimation.
in comparison bratnextgen identifies recombination by finding patterns of shared snps and so should perform better on this diverse dataset, but did require the analysis of the clades separately.
search for the cause of the population partition in h. parasuis
evidence of different restriction modification systems modulating homologous recombination has been associated with clades of neisseria meningitidis but evidence of recombination between the clades was still found in this bacterium  <cit> .
several methods were used to try to elucidate the differences between the baps populations and the two clades of h. parasuis using genome wide association methods , applied to both the core and the accessory genomes.
as we have seen some separation of serovars between these clades , it is possible that differences in the expression patterns of the capsule between serovars may present a recombination barrier  <cit> .
however, recombination within the clades is still occurring between isolates of different serovars.
the differences between the clades based on the accessory genome and the phylogeny  was visualised by creating heat-maps of the shared accessory genes between isolates.
the heat-map is ordered based on the phylogeny which was built using the snps from the core genome , and a dendrogram built based on the similarity in the accessory genome or shared genes .
clear separation can be seen in the accessory genome between the two clades, shown in figure  <dig>  with only an additional 500– <dig> accessory genes shared between the two clades aside from the  <dig> genes from the core genome.figure  <dig>
heat-map of the shared accessory genes between strains.
the plot is ordered by the phylogeny based on the snps within the core genome .
while the columns are ordered based on the similarity in the presence and absence pattern of accessory genes between isolates, which is represented by a dendrogram along the top of the heat-map.
a clear separation can be seen between the clades, and it appears that both the phylogeny and the dendrogram split the population into the two clades, suggesting little recombination occurs between the two clades, but there is recombination within them both.
additionally, the separation of two clades based on the accessory genome strongly suggests that horizontal transfer occurs in both clades, but that transfer of genetic material is restricted between these groups.
this plot allowed us to identify several areas of high similarity between strains, for example isolates hs <dig>  hs <dig>  hs <dig>  hs <dig>  hs <dig>  hs <dig> and hs <dig> share five intact phage regions that were previously identified in aeromonas , burkholderia , mannheimia , salmonella  and shigella  based on their hits from the phast database.
furthermore, a plasmid with similarity to one isolated from haemophilus influenzae  was identified in nine isolates using iceberg.
in addition, we were able to identify an integrated antibiotic resistance gene element with greater than 99% identity to one identified in vibrio cholera , in the two strains possessing multiple antibiotic resistance genes .
a blast analysis of all known h. parasuis genes involved in recombination, competence and restriction modification genes was performed using the isolate collection, but the only difference found was an overrepresentation of transposases in clade  <dig> compared to clade  <dig>  experimental data on the natural competence of isolates in this collection are sparse.
we have tested a small number of isolates to date, including strains  <dig> , hs <dig>  and hs <dig> , which have been classed as naturally transformable, while  <dig> is not transformable .
the presence of naturally competent isolates in both clades suggests that the partition between the two clades is not due to differences in competency genes.
therefore, to investigate this separation between the clades, we performed a gwas using discriminant analysis of principal components  on the baps populations.
this allowed us to identify genotypes  associated with population structure  within our isolate collection.
we also used generalised linear modelling  to test if these genotypes were significantly different between the baps populations , whether snp-based or based on gene presence/absence.
we varied the quantities of pca eigen-values  included in the dapc as a cross-validation technique to allow for comparison of the candidates that correlated with population structure.
the dapc analysis of the core genome showed a clear separation between the baps populations , and this is also evident between the two clades .
similarly, separation of the baps populations and clades can be seen based on the accessory genome  although to a lesser extent than from the core genome.
due to the large separation between the baps populations and the clades, a large number of genes and snps were identified from this analysis.
even with the top 10% of the loadings from the dapc ; the numbers of genotypes identified from the dapc cross-validation was greater than  <dig> for both snps and the accessory genome.
the shortlist of  <dig> snps identified from the dapc of the core genome is shown in additional file 10: table s <dig>  however, only  <dig> snps were identified as significantly different between baps populations or the two clades and they did not result in a predicted amino acid change in the encoded protein .
two snps were identified within the holliday junction atp-dependent dna helicase  and resulted in amino acid changes.
it is interesting to speculate that these snps may explain the difference in recombination rates between the two clades.
the shortlist of  <dig> genes identified from the dapc of the accessory genome is shown in additional file 11: table s <dig>  with only four genes identified as significantly different between baps populations .
of the genes in the shortlist, there are several genes encoding outer membrane proteins and potential capsule polysaccharide proteins, two transcriptional regulators, and genes of unknown function  that may potentially play a role in the uptake or exchange of dna.
pfam analysis  <cit>  of the proteins of unknown function  identified domains involved in transporter-associated domains, a helix-turn-helix domain and a lipoprotein domain, and three domains of unknown function, which are not obviously associated with the separation of the two clades.
based on the heat-map  and the dapc plots , we suggest that there are a greater number of genes that are contributing to this separation but they are not found exclusively in particular baps populations or clades.
this analysis has helped to understand some of the differences between the groups, and the snps within the holliday junction helicase are worth investigating further.
the pronounced difference between the two clades of h. parasuis warrants further investigation to help us to understand the separation of the two lineages in h. parasuis at both a genomic and functional level and their potential involvement in the barrier to recombination between these two clades.
genomic differences between serovars
historically, serotyping has been the main basis of subcategorising h. parasuis, and the currently available vaccines are serovar-specific.
although surveillance is not carried out routinely in many countries, it appears that the epidemiological profile of the serovars has been fairly stable throughout the last 30 years  <cit> .
multiple sets of reference strains of h. parasuis are currently in use around the world, and so we have sequenced these strains and confirmed that those with the same name are identical, with the exception of d <dig> .
the two strains labelled d <dig> shared only  <dig> genes of  <dig> genes  and differ by  <dig> snps in the core genome, which places them within the same clade but they belong to different baps populations.
reference strains that differed between sets were found to possess the same capsule locus but were otherwise distantly related based on the core and accessory genome.
we have found a high association between the capsule locus types and the serotyped isolate collection and propose that the capsule locus is the main determinant of serovar as has been previously assumed  <cit>  .
this association was also evident from capsule locus synteny  when the organisation of the capsule locus of the reference strains was compared by blastn to the remainder of the isolate collection  <cit> .
we have also examined the serotyped isolate collection using dapc and glm to try to identify genes involved in serovar designation, this being detailed in additional file  <dig> and additional file 13: figure s <dig>  briefly, based on the prevalent serovars in the isolate collection we were able to identify twenty-three potential serovar-specific genes, five of these candidates were previously identified capsule genes, three phage genes, a transposase, a filamentous haemagglutinin and  <dig> have unknown functions.
this is the first time that a gwas method has been used to try understand an existing typing scheme.
this analysis was based on a relatively small subset of strains when considering the number of isolates of each individual serovar, and so with a greater number of the less prevalent serovars, these methods have a higher likelihood of finding further serovar-specific gene markers.
this would allow for the refinement of the typing scheme to a molecular based test, which could have higher turn-around time, better sensitivity and specificity, and maximal reproducibility between laboratories through not depending on antibody reagents and the technical prowess of the laboratory worker.
genomic differences in the pan-genome by disease association
putative virulence factors are mostly found in the core genome
the majority of studies proposing virulence factors for this bacterium have focused on a subset of the most prevalent serovars or on a small number of isolates, and so we investigated the distribution of these suggested virulence factors within our more diverse isolate collection.
a survey of the literature produced a list of  <dig> putative virulence factors that were retrievable from ncbi as either a nucleotide or protein sequence  .
this included the vtaa genes, which are some of the best characterized virulence markers for h. parasuis identified to date  <cit> , but which were found at low frequency in this isolate collection.
this may be due to the repetitive nature of the nucleotide sequences of trimeric auto-transporters, which means that these may not have assembled well in our draft genome data.
instead we used the translocator and leader sequences of these genes and identified the group  <dig>   <dig> and  <dig> translocator domains in all isolates.
of these proposed virulence factors,  <dig> were present at greater than 80% identity in all of our isolates and so were present in the core genome.
there is the possibility that variation within these genes may contribute to differences in virulence but this analysis was beyond the scope of this study.
the distribution of the putative virulence factors from the accessory genome was examined with respect to disease association  and disease phenotypes  by the calculation of their relative proportions.
five virulence factors were over-represented in the non-clinical isolates  and under-represented in the clinical isolates.
conversely,  <dig> virulence factors were over-represented in the clinical isolates .
when considering the disease phenotypes, no virulence factors were over-represented in only the respiratory population and only hsdr and f357_gp <dig> were over-represented in systemic isolates.
we also looked into the differences in genome size and number of genes between the disease associated isolates, as a proxy for evidence of reductive evolution  <cit> .
no evidence of reductive evolution was found based on our isolate collection; details can be found in additional file  <dig> and additional file  <dig>
a high number of genes from the accessory genome were linked to differences in clinical phenotype
while we were able to find some associations between existing putative virulence factors for h. parasuis and clinical phenotype, there is also a large amount of data from the pan-genome that could yield new important snps or genes that are associated with pathogenicity.
therefore, we used the methods of dapc and glm to assess the relationship between the pan-genome and virulence based on the clinical metadata.
from this analysis we have identified both snps and genes that allow us to separate the isolate collection into disease categories.
all of the candidates were added into the glm with clinical disease as the factor in the model, and those that were significantly different between the clinical and non-clinical isolates  were considered as potential markers for virulence.
we also performed analyses of disease phenotypes for both the core and the accessory genome, but very little difference between the respiratory and systemic isolates was identified using dapc .
this occurred whether isolates from the lung with signs of systemic disease were classed as either respiratory or systemic strains.
this suggests that the difference between respiratory or systemic outcome of disease may be dependent on additional host factors such as welfare and host immune response rather than genetic differences between the isolates or particular tissue tropisms.
the discriminant function of the disease association is shown in figures 5a and b for the core and accessory genomes respectively .
this shows separation between the clinical and non-clinical isolates to a similar extent for both the core and the accessory genome.
some overlap was seen between these categories , but this is not unexpected as some isolates may be capable of causing disease but were not doing so at the time of isolation.
we also considered the influence of geography and the separation of the strains into populations with different recombination rates  on the pan-genome and saw that both showed separation in the discriminant functions of the core and accessory genomes .figure  <dig>
discriminant analysis of principal components applied to clinical strains of
h. parasuis
by disease categories.
80% eigenvalues were retained for the pca and all eigenvalues were retained for the discriminant analysis.
plots a and b show the first axis of the discriminant function while c, d, e and f show the first two axes.
separation along the axes suggests that genetic differences are present between the phenotypic groups of clinical and non-clinical isolates that are being compared; however the presence of overlap shows that some strains are intermediates.
plots c and d show that geographic origin, do not show much separation of the isolates by geography, those that have separated are only represented by a couple of isolates.
on the other hand, the discriminant function based on the baps populations shows a lot of separation and so the population structure does have an influence on these isolates genetic content.
when we consider the candidate markers for virulence from the dapc and the glm, a greater number of candidates were significantly different between the clinical and non-clinical isolates from the glm of the accessory genome  over the core .
of the snps identified as being statistically significant, seven resulted in predicted amino acid changes.
however, only four predicted the change of amino acid type, which may influence the structure of the protein.
for example, two snps identified in a long chain fatty acid coa ligase  and in a transcriptional regulatory protein , which would result in an a/s amino acid change and a change from hydrophobic to polar.
a snp in a fatty acid oxidation complex subunit alpha  was predicted changing the amino acid from hydrophobic to polar and finally a aromatic to hydrophobic change at holo-acyl carrier protein synthase .
a gene involved in fatty acid degradation  was previously suggested as having a role in virulence in h. parasuis <cit> .
furthermore, the association of fatty acid metabolism and virulence has been previously identified in mycobacterium tuberculosis <cit> .
of the genes yielded from the accessory genome analysis nine were putative virulence factors from the previous analysis of published virulence factors .
from the dapc pathogenicity candidates we identified genes encoding four adhesins, abc transporters, an antirepressor, an aspartokinase, two cytolethal distending toxins, enoyl-coa hydratase/carnithine racemase-like protein, ferric uptake proteins, glycerol uptake facilitator protein, a hemolysin transporter, a methylase, an n-acetylmuramoyl-l-alanine amidase, a nucleoside symporter, an outer membrane protein, four phage proteins, an rna helicase, a putative serine protease, a tonb dependent receptor, another toxin , a transcriptional activator, two transposases, a xylose import atp-binding protein and nineteen proteins of unknown function.
the genes of unknown function were run through the pfam database  to try to elucidate any further details about their function, but only three had any hits to known domains, one with a phage domain , one with an abc transporter and membrane domains  and another with a yada adhesin domain .table  <dig>
significantly different genes identified from h. parasuis between clinical and non-clinical isolates using dapc
h. parasuis
aspartokinase lysc
atpdependent rna helicase hrpa
modification methylase hpaii
putative outer membrane usher protein htre
transcriptional activator prtr
xylose import atpbinding protein xylg
all genes were identified as significant using generalised linear modelling of different dapc iterations using between 60-90% of the principal components .
the phenotypic data used in the model was whether isolates were clinical or non-clinical compared to the presence or absence of the genes.
these genes are in the same families as other common virulence determinants in other better-studied pathogens.
for example, adhesins are commonly regarded as virulence factors as they allow adhesion to the host cell, and depending on their specificity may be required for the progression of disease in particular tissues  <cit> .
while cytolethal distending toxins have been previously studied in h. parasuis <cit> , there appear to be multiple, divergent copies within the isolate collection.
variation in these cdt genes has previously been found between different strains of e. coli at as little at 39% amino acid identity, and so the presence of multiple copies of this family of toxins may impact the virulence of isolates of h. parasuis <cit> .
in addition, glycerol uptake facilitators have been implicated in virulence and growth in mycoplasma pneumoniae <cit> ; n-acetylmuramoyl-l-alanine amidases play a role in the conversion of the streptococcus mitis autolysin  <cit> ; serine proteases are regarded as virulence factors in enterohemorrhagic e. coli <cit>  and streptococcus pneumoniae <cit>  and the tonb-dependent receptors are involved in iron acquisition, which is important in pathogens such as neisseria meningitidis <cit>  and vibrio cholera <cit> .
in summary, the dapc and glm analyses of the isolate collection have yielded a short-list of new putative virulence factors for further investigation in this pathogen and their use as epidemiological markers will be considered further.
CONCLUSIONS
we have investigated the composition of the pan-genome of the major pig pathogen h. parasuis with respect to detailed metadata in order to investigate associations between genotype and phenotype.
we have found that this is a diverse bacterial species, comprised of two main lineages, exhibiting very high rates of recombination, which appears to be restricted within one or other of the two clades.
using dapc as the gwas method, we were able to identify snps within a holliday junction atp-dependent dna helicase  that may be involved in the partitioning of this species into two clades.
we have found genetic differences with regard to serovar based on the capsule loci and further potential serovar-specific genes from the accessory genome, which is the first time that a gwas method has been used to try to refine an existing typing method.
we have also identified genetic differences in the clinical and non-clinical isolates from  <dig> genes, nine of which were also identified by dapc and glm analyses.
the discriminant power of the accessory genome over the core genome provides a better list of candidate virulence genes that will be invaluable for investigating the pathogenesis of h. parasuis and how it causes glässer’s disease.
the accessory genome is therefore an important source for understanding the pathogenesis of diverse species of bacteria, particularly in combination with detailed metadata.
these methods are applicable to any pathogen where detailed clinical metadata can be obtained and analysed in combination with next generation sequencing, which could identify genotypes  that may be useful in the generation of new molecular diagnostic and vaccine candidates that can aid surveillance and prevention of disease.

 haemophilus parasuis is an important pathogen of pigs but is also a frequent commensal of their upper respiratory tracts. it causes polyserositis, meningitis and septicaemia as well as pneumonia and pleurisy. historically, h. parasuis caused disease in recently weaned piglets, but with the intensification of production it now affects a wide age-range. h. parasuis is the leading cause of mortality, alongside the prrs virus, in nursery herds in the usa, and is the third-most important bacterial pathogen affecting finisher herds. h. parasuis also contributes to multi-factorial porcine respiratory disease complex, the leading cause of mortality in grower-finisher pigs in the usa. h. parasuis is endemic in all pig farming countries and has a major economic impact on the global pig industry, with costs including those for production and stock losses, carcass disposal, and vaccination, as well as requiring antibiotic therapy. the study of h. parasuis to date has classified the bacterium, originally using a gel immuno-diffusion assay and more commonly now an indirect haemagglutination assay, into different serovars with non-typeable isolates also frequently found. the designation of serovar in h. parasuis is thought to be predominantly due to the presence of a particular polysaccharide capsule which is specific to each serovar. some serovars, especially and cause disease more commonly than others. however, no absolute relationship between virulence and serovar has been found. the onset of disease has been linked to many factors, including physiological stress and immune status, but the current study focuses on understanding the bacterial genetic factors that might correlate with the establishment of disease. several experiments have tried to associate the bacterium’s genotype with phenotype, but with limited success. experimental evidence to support definitive virulence factors in h. parasuis is sparse and most studies have been limited to small numbers of isolates or serovars. it is likely that the relationship between genotype and virulence phenotype is complex. the first h. parasuis genome to be sequenced was in, with several further sequences released since then. these have been used to design several experiments, including searches for novel immunogenic proteins and virulence factors, but these few strains may not be reflective of the wider diversity of isolates that exist in pig populations. previous studies show that the pan-genome of a bacterial species can be vast. recombination and lateral transfer are involved in creating the pan-genome, which can be resolved into the core and accessory genomes. the boundaries of the core genome can be extrapolated from a small number of genomes, but it is preferable to analyse a large number of isolates for a more comprehensive definition. an extensive amount of variation can often be found within the accessory genome of a species and it is important to investigate this source of variation to understand the adaptive potential of a given species. reductions in cost have allowed the whole genome sequencing of many strains of the same species, and this can provide the necessary statistical power required to find associations between genotypes and phenotypes such as virulence. this variation is now being explored in many pathogens including campylobacter jejuni, escherichia coli, haemophilus influenzae, staphylococcus aureus and streptococcus pneumoniae. such studies seek associations between phenotypes and genotypes that would inform surveillance efforts, improve diagnostics, and aid understanding of the pathogen, potentially yielding new treatments of the associated diseases. in order to disentangle the genetic and environmental bases of glässer’s disease, we undertook a pan-genome analysis of isolates of h. parasuis with high quality clinical metadata, covering all serovars and isolates from disease- and non-disease-causing backgrounds. we have investigated the key features of the pan-genome and assessed the relationship between the genetics of this bacterium and collated metadata, using statistical analyses of single nucleotide polymorphisms in the core genome, and the presence of genes in the accessory genome. we compared our analyses to the distribution of previously published virulence factors of h. parasuis within our collection of isolates. we found that the composition of the accessory genome was a significant factor in determining whether isolates were likely to cause disease. our findings emphasise the importance of describing the pan-genome of a bacterial species in order to understand pathogenesis, and to underpin a new generation of vaccine and diagnostic candidates for endemic bacterial pathogens.  epidemiology of the isolate collection the isolates in the collection originate from the uk, argentina, belgium, china, denmark, germany, greece, italy, japan, spain, sweden, switzerland and the usa, with nine isolates of unknown origin. our analysis included the reference strains of the serovars and the published genomes of the us strain and the chinese strains sh and zj. the dates of collection ranged between 1991– for european isolates and 1955– for the reference strains. for the uk isolates, we collected clinical metadata including veterinary diagnostic information from necropsy, the tissue of isolation and details of any other signs of disease in the pig. clinical metadata for the other strains included tissue of origin, clinical signs and, for a subset of strains, experimental reproduction of disease. we used this information to categorise the isolates into clinical and non-clinical categories. a non-clinical isolate was an isolate identified in the upper respiratory tract in the absence of any signs of respiratory or systemic disease typical of infection with h. parasuis. a clinical isolate was an isolate identified from an animal displaying signs of h. parasuis-related disease and isolated from a respiratory or a systemic site. we divided the clinical isolates further into respiratory and systemic, to reflect disease severity. respiratory isolates were obtained from the lung without signs of systemic disease elsewhere in the animal. systemic isolates were obtained from a systemic site e.g. brain, joint or blood which are sites typically affected in glässer’s disease. isolates that were taken from the lung of an animal that also showed signs of systemic disease were first tested in the statistical analyses as respiratory isolates and then as systemic. this was to allow for the possibility of these isolates being from a co-infection in the same animal or where the same isolate caused disease elsewhere in the same animal. one hundred and seventeen of the isolates had been previously serotyped and these included examples of all fifteen known serovars. for details of the isolate collection including disease categories and serotyping results see additional file 1: table s analysis of the distribution of isolates by serovar and disease association showed that the most prevalent serovar in this collection was serovar followed by non-typeables and serovars and which fits with global prevalence data. our data support the previous limited studies that show that certain serovars, particularly and were more strongly linked to disease. while serovar is still important for the formulation of vaccine strategy against this bacterium, its relationship with virulence has long been discussed but never proven. the fact that there are non-clinical isolates within the same serovars highlights the limitation and uncertainty of serotyping as a technique for predicting disease potential. we expect that there are isolates amongst the non-clinical collection that will have the potential to cause clinical disease given a different host or environment, and also that some of the clinical isolates caused disease due to other environmental factors or co-infections. the main reason for this study was to investigate the balance between the host, environment and genotypes for example serovar determining factors, predicted virulence genes and to discover new gene associations with phenotypes such as virulence. we propose that molecular epidemiological markers that relate to severity of disease would be more useful for h. parasuis surveillance now that rapid genome sequencing of isolates is possible. to achieve this goal we undertook pan-genome and statistical analyses with the two-fold aim of identifying genes that correlate with virulence and so could be worth investigating as markers of disease-causing potential, and with serovar to allow the design of a molecular-based serotyping technique.figure h. parasuis isolate collection displayed by serovar and disease. each serovar includes those found with cross-reactions. nt denotes'non-typeable’ by serotyping either due to no reaction or due to more than three reactions to the serotyping antisera. strains without serotyping data are excluded. the pan-genome of h. parasuis we use the term pan-genome to refer to the set of all of the genes predicted from every isolate in the collection. the pan-genome comprised homology groups encoding genes from the orthomcl analysis of the isolates. of these, genes were found within all isolates and so were designated as the core genome. this core genome is approximately  mb of the  mb average genome size with snps identified. the accessory genome was made up of homology groups that varied between isolates. of these homology groups were classed as pseudogenes, and % of these homology groups were identified in less than 1% of the isolates in the collection. for details of the cog functional classification of the pan- genome see additional file and additional file 3: figure s synteny of the pan-genome the syntenic pan-genome of h. parasuis is the pan-genome reordered to reflect the chromosomal location of each gene based on the sh strain. genes were inserted from contigs in descending size order, based on their order along each contig. as an iterative process not all genes could be assigned to a location in the syntenic pan-genome. therefore the syntenic pan-genome of h. parasuis represents the genes out of the total genes where the position could be identified in the pan-genome order, and is shown in additional file 4: figure s several regions of variation can be seen within additional file 4: figure s with the location of many phage genes, highly variable regions and antibiotic resistance genes denoted by colouring of the heat-map. the capsule genes did not cluster in one location in the syntenic pan-genome, which may be due to recombination in the genomes, and so the synteny of the capsule locus was established separately and can be seen in additional file 5: figure s several phages were identified within this collection including the mu phage, a siphovirus as well as other phage genes throughout the pan-genome. for details of genes identified within these variable regions see additional file briefly, within these regions we identified many phage genes, two conjugal transfer operons, adhesins, toxins, iron-sulphur binding proteins and many proteins of unknown function. based on the similarity of the g + c content to the average for this bacterium, many of the genes in region appear to be from h. parasuis genes that have yet to be characterised and submitted to ncbi. regions and contained some genes of relatively high g + c content but no hits were found when blast analysis was performed for these encoded proteins using the ncbi non-redundant protein database. only region contained multiple proteins with hits from nr to other bacterial species including actinobacillus, gallibacterium, klebsiella and mannhaimia. the average g + c content of the pan-genome was % g + c content, ranging between - % g + c content for all isolates, which is discussed further in additional file and additional file 6: figure s population structure of h. parasuis identification of two divergent clades of h. parasuis analyses of the population structure of this isolate collection was performed using bayesian analysis of population structure on the core genome. five baps populations were identified from these analyses, these represented groups within the overall population that share similar snps and indicate a common ancestry with similar amounts of admixture. the baps plot shown in additional file 7: figure s demonstrates the within-cluster and between-cluster variation for the five baps populations based on the core genome. while the partitioning of population structure in other bacteria has been reported before, further analysis of the causes of the partitioning usually leads to an explanation dependent on ecological separation e.g. due to geographical separation or restriction to different hosts or tissues. h. parasuis has only been identified in pigs and wild boar, and the route of infection is via the nose and the upper respiratory tract and so are all able to inhabit the same niche. furthermore, the global nature of the pig industry means that import and export of stock occurs worldwide and so we did not expect that geographical separation would have much influence on the population structure. we began our investigation into the population structure by checking these potential ecological barriers by comparing our clinical metadata to the phylogenetic tree of the h. parasuis core genome built from the concatenated core genes. the topology of the phylogeny shows the segregation of the tree into two clades, which fits with the separation of the baps populations with baps populations and found in clade and groups and in clade we also investigated the correlation between the baps populations and the clinical metadata. when both the core genome and baps populations were compared to the metadata, some separation of the serovars can be seen between the two clades, with serovars and identified predominantly in clade serovar is the only serovar to be found in equal proportions in both the clades. serovar is mostly found in clade for the remainder of the serovars, the number of isolates was too low to be able to restrict them to one clade. in comparison, no influence of disease association on population structure was observed, as the clinical isolates were found throughout the tree and the baps populations. this is contradictory to the population grouping predicted by the multi-locus sequence typing scheme, which showed six main subgroups on the mlst phylogeny; one predominantly clinical, one mostly non-clinical and four mixed clades. analysis of the individual baps populations in figure  showed that baps populations and contained very few non-clinical isolates. some geographic structure can be seen in figure  with a greater proportion of the uk isolates in clade which is likely to be due to the large proportion of uk isolates in the collection. with the global nature of the pig market it is not unexpected that isolates from different geographic locations or serovars would be found in each of the baps populations. indeed, populations and are very mixed; indicating only limited separation by a geographic barrier. the baps populations are also very mixed by tissue of isolation and so no tissue tropism can be seen based on the population structure. of the categories of metadata, only serovar shows some association with the baps populations and the differences in recombination. therefore we conducted a second level baps analysis which split the current five baps populations into more defined populations; however the second level baps analysis still did not explain the population structure.figure core genome neighbor-joining tree of h. parasuis. trees are overlaid with the populations from the bayesian analysis of population structure, which represent isolates with similar rates of homologous recombination. further metadata including disease association, serovar and country of origin are also shown. baps populations explain the separation of the isolates on the tree into two main clades. the baps populations defined using bayesian analysis of population structure represent isolates with similar rates of homologous recombination. the baps populations have been compared to serovar, disease association and country of origin. some similarity in the baps populations and can be seen based on serovar distribution. very little difference can be seen between baps population when disease association is considered. uk isolates are overrepresented in baps populations and while and are of mixed geographic origin. high rates of recombination were identified in h. parasuis analysis of the regions undergoing recombination in the core genome was performed using bayesian recombination tracker on the two separate clades present in figure  we decided to perform separate analyses because the long-branch separating clades and would reduce the statistical power to detect recombination in either clade in a joint analysis. the recombination analysis was also performed using gubbins, which was performed on the alignment for all isolates. we were able to calculate the recombination rate relative to the mutation rate or homologous recombination rate for the individual baps populations using the output from gubbins; baps had the highest r/m of followed by baps, baps, baps and finally baps. however, each baps population contained isolates with a range of recombination rates, with multiple isolates with a r/m of when the second level baps populations were considered, some clusters of isolates with r/m of were identified but the majority of these isolates were still separated between the groups identified in the second level analysis. even with these isolates with low estimates of recombination rates, the r/m for each baps population is very high compared to published recombination rates in bacteria and is far higher than previously estimated for h. parasuis. one previous study calculated the r/m for h. parasuis based on the genes involved in the mlst scheme, and identified higher rates for the clinical isolates compared with the non-clinical isolates. our estimates for the rates of homologous recombination in h. parasuis far exceed the estimates based on the mlst housekeeping genes. this may be due to purifying selection that is likely to be affecting housekeeping genes and so the estimate based on the core genome, which contains far more genetic information, should give a more informative estimate of recombination across the genome. the values of r/m of the two clades of h. parasuis estimated using both recombination analysis methods were within an order of magnitude. the r/m for the first clade was and, while the r/m for the second clade was estimated at and, using gubbins and bratnextgen respectively. while both methods yielded comparable estimates of homologous recombination rates, the diversity in this h. parasuis population is greater than is usually included in the analysis using gubbins and so care should be taken when interpreting the results. gubbins was designed to identify recombinations from distant sources within very clonal populations and so these r/ms may even be an under-estimation. in comparison bratnextgen identifies recombination by finding patterns of shared snps and so should perform better on this diverse dataset, but did require the analysis of the clades separately. search for the cause of the population partition in h. parasuis evidence of different restriction modification systems modulating homologous recombination has been associated with clades of neisseria meningitidis but evidence of recombination between the clades was still found in this bacterium. several methods were used to try to elucidate the differences between the baps populations and the two clades of h. parasuis using genome wide association methods, applied to both the core and the accessory genomes. as we have seen some separation of serovars between these clades, it is possible that differences in the expression patterns of the capsule between serovars may present a recombination barrier. however, recombination within the clades is still occurring between isolates of different serovars. the differences between the clades based on the accessory genome and the phylogeny was visualised by creating heat-maps of the shared accessory genes between isolates. the heat-map is ordered based on the phylogeny which was built using the snps from the core genome, and a dendrogram built based on the similarity in the accessory genome or shared genes. clear separation can be seen in the accessory genome between the two clades, shown in figure  with only an additional 500– accessory genes shared between the two clades aside from the genes from the core genome.figure heat-map of the shared accessory genes between strains. the plot is ordered by the phylogeny based on the snps within the core genome. while the columns are ordered based on the similarity in the presence and absence pattern of accessory genes between isolates, which is represented by a dendrogram along the top of the heat-map. a clear separation can be seen between the clades, and it appears that both the phylogeny and the dendrogram split the population into the two clades, suggesting little recombination occurs between the two clades, but there is recombination within them both. additionally, the separation of two clades based on the accessory genome strongly suggests that horizontal transfer occurs in both clades, but that transfer of genetic material is restricted between these groups. this plot allowed us to identify several areas of high similarity between strains, for example isolates hs hs hs hs hs hs and hs share five intact phage regions that were previously identified in aeromonas, burkholderia, mannheimia, salmonella and shigella based on their hits from the phast database. furthermore, a plasmid with similarity to one isolated from haemophilus influenzae was identified in nine isolates using iceberg. in addition, we were able to identify an integrated antibiotic resistance gene element with greater than 99% identity to one identified in vibrio cholera, in the two strains possessing multiple antibiotic resistance genes. a blast analysis of all known h. parasuis genes involved in recombination, competence and restriction modification genes was performed using the isolate collection, but the only difference found was an overrepresentation of transposases in clade compared to clade experimental data on the natural competence of isolates in this collection are sparse. we have tested a small number of isolates to date, including strains, hs and hs, which have been classed as naturally transformable, while is not transformable. the presence of naturally competent isolates in both clades suggests that the partition between the two clades is not due to differences in competency genes. therefore, to investigate this separation between the clades, we performed a gwas using discriminant analysis of principal components on the baps populations. this allowed us to identify genotypes associated with population structure within our isolate collection. we also used generalised linear modelling to test if these genotypes were significantly different between the baps populations, whether snp-based or based on gene presence/absence. we varied the quantities of pca eigen-values included in the dapc as a cross-validation technique to allow for comparison of the candidates that correlated with population structure. the dapc analysis of the core genome showed a clear separation between the baps populations, and this is also evident between the two clades. similarly, separation of the baps populations and clades can be seen based on the accessory genome although to a lesser extent than from the core genome. due to the large separation between the baps populations and the clades, a large number of genes and snps were identified from this analysis. even with the top 10% of the loadings from the dapc; the numbers of genotypes identified from the dapc cross-validation was greater than for both snps and the accessory genome. the shortlist of snps identified from the dapc of the core genome is shown in additional file 10: table s however, only snps were identified as significantly different between baps populations or the two clades and they did not result in a predicted amino acid change in the encoded protein. two snps were identified within the holliday junction atp-dependent dna helicase and resulted in amino acid changes. it is interesting to speculate that these snps may explain the difference in recombination rates between the two clades. the shortlist of genes identified from the dapc of the accessory genome is shown in additional file 11: table s with only four genes identified as significantly different between baps populations. of the genes in the shortlist, there are several genes encoding outer membrane proteins and potential capsule polysaccharide proteins, two transcriptional regulators, and genes of unknown function that may potentially play a role in the uptake or exchange of dna. pfam analysis of the proteins of unknown function identified domains involved in transporter-associated domains, a helix-turn-helix domain and a lipoprotein domain, and three domains of unknown function, which are not obviously associated with the separation of the two clades. based on the heat-map and the dapc plots, we suggest that there are a greater number of genes that are contributing to this separation but they are not found exclusively in particular baps populations or clades. this analysis has helped to understand some of the differences between the groups, and the snps within the holliday junction helicase are worth investigating further. the pronounced difference between the two clades of h. parasuis warrants further investigation to help us to understand the separation of the two lineages in h. parasuis at both a genomic and functional level and their potential involvement in the barrier to recombination between these two clades. genomic differences between serovars historically, serotyping has been the main basis of subcategorising h. parasuis, and the currently available vaccines are serovar-specific. although surveillance is not carried out routinely in many countries, it appears that the epidemiological profile of the serovars has been fairly stable throughout the last 30 years. multiple sets of reference strains of h. parasuis are currently in use around the world, and so we have sequenced these strains and confirmed that those with the same name are identical, with the exception of d. the two strains labelled d shared only genes of genes and differ by snps in the core genome, which places them within the same clade but they belong to different baps populations. reference strains that differed between sets were found to possess the same capsule locus but were otherwise distantly related based on the core and accessory genome. we have found a high association between the capsule locus types and the serotyped isolate collection and propose that the capsule locus is the main determinant of serovar as has been previously assumed. this association was also evident from capsule locus synteny when the organisation of the capsule locus of the reference strains was compared by blastn to the remainder of the isolate collection. we have also examined the serotyped isolate collection using dapc and glm to try to identify genes involved in serovar designation, this being detailed in additional file and additional file 13: figure s briefly, based on the prevalent serovars in the isolate collection we were able to identify twenty-three potential serovar-specific genes, five of these candidates were previously identified capsule genes, three phage genes, a transposase, a filamentous haemagglutinin and have unknown functions. this is the first time that a gwas method has been used to try understand an existing typing scheme. this analysis was based on a relatively small subset of strains when considering the number of isolates of each individual serovar, and so with a greater number of the less prevalent serovars, these methods have a higher likelihood of finding further serovar-specific gene markers. this would allow for the refinement of the typing scheme to a molecular based test, which could have higher turn-around time, better sensitivity and specificity, and maximal reproducibility between laboratories through not depending on antibody reagents and the technical prowess of the laboratory worker. genomic differences in the pan-genome by disease association putative virulence factors are mostly found in the core genome the majority of studies proposing virulence factors for this bacterium have focused on a subset of the most prevalent serovars or on a small number of isolates, and so we investigated the distribution of these suggested virulence factors within our more diverse isolate collection. a survey of the literature produced a list of putative virulence factors that were retrievable from ncbi as either a nucleotide or protein sequence. this included the vtaa genes, which are some of the best characterized virulence markers for h. parasuis identified to date, but which were found at low frequency in this isolate collection. this may be due to the repetitive nature of the nucleotide sequences of trimeric auto-transporters, which means that these may not have assembled well in our draft genome data. instead we used the translocator and leader sequences of these genes and identified the group and translocator domains in all isolates. of these proposed virulence factors, were present at greater than 80% identity in all of our isolates and so were present in the core genome. there is the possibility that variation within these genes may contribute to differences in virulence but this analysis was beyond the scope of this study. the distribution of the putative virulence factors from the accessory genome was examined with respect to disease association and disease phenotypes by the calculation of their relative proportions. five virulence factors were over-represented in the non-clinical isolates and under-represented in the clinical isolates. conversely, virulence factors were over-represented in the clinical isolates. when considering the disease phenotypes, no virulence factors were over-represented in only the respiratory population and only hsdr and f357_gp were over-represented in systemic isolates. we also looked into the differences in genome size and number of genes between the disease associated isolates, as a proxy for evidence of reductive evolution. no evidence of reductive evolution was found based on our isolate collection; details can be found in additional file and additional file a high number of genes from the accessory genome were linked to differences in clinical phenotype while we were able to find some associations between existing putative virulence factors for h. parasuis and clinical phenotype, there is also a large amount of data from the pan-genome that could yield new important snps or genes that are associated with pathogenicity. therefore, we used the methods of dapc and glm to assess the relationship between the pan-genome and virulence based on the clinical metadata. from this analysis we have identified both snps and genes that allow us to separate the isolate collection into disease categories. all of the candidates were added into the glm with clinical disease as the factor in the model, and those that were significantly different between the clinical and non-clinical isolates were considered as potential markers for virulence. we also performed analyses of disease phenotypes for both the core and the accessory genome, but very little difference between the respiratory and systemic isolates was identified using dapc. this occurred whether isolates from the lung with signs of systemic disease were classed as either respiratory or systemic strains. this suggests that the difference between respiratory or systemic outcome of disease may be dependent on additional host factors such as welfare and host immune response rather than genetic differences between the isolates or particular tissue tropisms. the discriminant function of the disease association is shown in figures 5a and b for the core and accessory genomes respectively. this shows separation between the clinical and non-clinical isolates to a similar extent for both the core and the accessory genome. some overlap was seen between these categories, but this is not unexpected as some isolates may be capable of causing disease but were not doing so at the time of isolation. we also considered the influence of geography and the separation of the strains into populations with different recombination rates on the pan-genome and saw that both showed separation in the discriminant functions of the core and accessory genomes.figure discriminant analysis of principal components applied to clinical strains of h. parasuis by disease categories. 80% eigenvalues were retained for the pca and all eigenvalues were retained for the discriminant analysis. plots a and b show the first axis of the discriminant function while c, d, e and f show the first two axes. separation along the axes suggests that genetic differences are present between the phenotypic groups of clinical and non-clinical isolates that are being compared; however the presence of overlap shows that some strains are intermediates. plots c and d show that geographic origin, do not show much separation of the isolates by geography, those that have separated are only represented by a couple of isolates. on the other hand, the discriminant function based on the baps populations shows a lot of separation and so the population structure does have an influence on these isolates genetic content. when we consider the candidate markers for virulence from the dapc and the glm, a greater number of candidates were significantly different between the clinical and non-clinical isolates from the glm of the accessory genome over the core. of the snps identified as being statistically significant, seven resulted in predicted amino acid changes. however, only four predicted the change of amino acid type, which may influence the structure of the protein. for example, two snps identified in a long chain fatty acid coa ligase and in a transcriptional regulatory protein, which would result in an a/s amino acid change and a change from hydrophobic to polar. a snp in a fatty acid oxidation complex subunit alpha was predicted changing the amino acid from hydrophobic to polar and finally a aromatic to hydrophobic change at holo-acyl carrier protein synthase. a gene involved in fatty acid degradation was previously suggested as having a role in virulence in h. parasuis. furthermore, the association of fatty acid metabolism and virulence has been previously identified in mycobacterium tuberculosis. of the genes yielded from the accessory genome analysis nine were putative virulence factors from the previous analysis of published virulence factors. from the dapc pathogenicity candidates we identified genes encoding four adhesins, abc transporters, an antirepressor, an aspartokinase, two cytolethal distending toxins, enoyl-coa hydratase/carnithine racemase-like protein, ferric uptake proteins, glycerol uptake facilitator protein, a hemolysin transporter, a methylase, an n-acetylmuramoyl-l-alanine amidase, a nucleoside symporter, an outer membrane protein, four phage proteins, an rna helicase, a putative serine protease, a tonb dependent receptor, another toxin, a transcriptional activator, two transposases, a xylose import atp-binding protein and nineteen proteins of unknown function. the genes of unknown function were run through the pfam database to try to elucidate any further details about their function, but only three had any hits to known domains, one with a phage domain, one with an abc transporter and membrane domains and another with a yada adhesin domain.table significantly different genes identified from h. parasuis between clinical and non-clinical isolates using dapc h. parasuis aspartokinase lysc atpdependent rna helicase hrpa modification methylase hpaii putative outer membrane usher protein htre transcriptional activator prtr xylose import atpbinding protein xylg all genes were identified as significant using generalised linear modelling of different dapc iterations using between 60-90% of the principal components. the phenotypic data used in the model was whether isolates were clinical or non-clinical compared to the presence or absence of the genes. these genes are in the same families as other common virulence determinants in other better-studied pathogens. for example, adhesins are commonly regarded as virulence factors as they allow adhesion to the host cell, and depending on their specificity may be required for the progression of disease in particular tissues. while cytolethal distending toxins have been previously studied in h. parasuis, there appear to be multiple, divergent copies within the isolate collection. variation in these cdt genes has previously been found between different strains of e. coli at as little at 39% amino acid identity, and so the presence of multiple copies of this family of toxins may impact the virulence of isolates of h. parasuis. in addition, glycerol uptake facilitators have been implicated in virulence and growth in mycoplasma pneumoniae; n-acetylmuramoyl-l-alanine amidases play a role in the conversion of the streptococcus mitis autolysin; serine proteases are regarded as virulence factors in enterohemorrhagic e. coli and streptococcus pneumoniae and the tonb-dependent receptors are involved in iron acquisition, which is important in pathogens such as neisseria meningitidis and vibrio cholera. in summary, the dapc and glm analyses of the isolate collection have yielded a short-list of new putative virulence factors for further investigation in this pathogen and their use as epidemiological markers will be considered further.  we have investigated the composition of the pan-genome of the major pig pathogen h. parasuis with respect to detailed metadata in order to investigate associations between genotype and phenotype. we have found that this is a diverse bacterial species, comprised of two main lineages, exhibiting very high rates of recombination, which appears to be restricted within one or other of the two clades. using dapc as the gwas method, we were able to identify snps within a holliday junction atp-dependent dna helicase that may be involved in the partitioning of this species into two clades. we have found genetic differences with regard to serovar based on the capsule loci and further potential serovar-specific genes from the accessory genome, which is the first time that a gwas method has been used to try to refine an existing typing method. we have also identified genetic differences in the clinical and non-clinical isolates from genes, nine of which were also identified by dapc and glm analyses. the discriminant power of the accessory genome over the core genome provides a better list of candidate virulence genes that will be invaluable for investigating the pathogenesis of h. parasuis and how it causes glässer’s disease. the accessory genome is therefore an important source for understanding the pathogenesis of diverse species of bacteria, particularly in combination with detailed metadata. these methods are applicable to any pathogen where detailed clinical metadata can be obtained and analysed in combination with next generation sequencing, which could identify genotypes that may be useful in the generation of new molecular diagnostic and vaccine candidates that can aid surveillance and prevention of disease. 
these genes will inform the generation of new molecular diagnostics and vaccines, and refinement of existing typing schemes and show the importance of the accessory genome of a diverse species when investigating the relationship between genotypes and phenotypes. we used genome-wide association methods including discriminant analysis of principal components and generalised linear modelling to look for genetic determinants of this population partition, serovar and pathogenicity. we were able to identify genes from the accessory genome that were significantly associated with phenotypes such as potential serovar specific genes including capsule genes, and putative virulence factors that were significantly different between the clinical and non-clinical isolates. we sequenced the genomes of isolates of h. parasuis and have used this to describe the pan-genome and to correlate this with clinical and carrier status, as well as with serotype. recombination and population structure analyses identified five groups with very high rates of recombination, separated into two clades of h. parasuis with no signs of recombination between them. haemophilus parasuis is the etiologic agent of glässer’s disease in pigs and causes devastating losses to the farming industry. genome wide association studypan-genomerecombinationdapcvirulence factorspopulation structurehaemophilus parasuisissue-copyright-statement© the author 2014
 haemophilus parasuis is the etiologic agent of glässer’s disease in pigs and causes devastating losses to the farming industry. whilst some hyper-virulent isolates have been described, the relationship between genetics and disease outcome has been only partially established. in particular, there is weak correlation between serovar and disease phenotype. we sequenced the genomes of isolates of h. parasuis and have used this to describe the pan-genome and to correlate this with clinical and carrier status, as well as with serotype.  recombination and population structure analyses identified five groups with very high rates of recombination, separated into two clades of h. parasuis with no signs of recombination between them. we used genome-wide association methods including discriminant analysis of principal components and generalised linear modelling to look for genetic determinants of this population partition, serovar and pathogenicity. we were able to identify genes from the accessory genome that were significantly associated with phenotypes such as potential serovar specific genes including capsule genes, and putative virulence factors that were significantly different between the clinical and non-clinical isolates. we also show that the presence of many previously suggested virulence factors is not an appropriate marker of virulence.  these genes will inform the generation of new molecular diagnostics and vaccines, and refinement of existing typing schemes and show the importance of the accessory genome of a diverse species when investigating the relationship between genotypes and phenotypes. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords genome wide association studypan-genomerecombinationdapcvirulence factorspopulation structurehaemophilus parasuisissue-copyright-statement© the author 2014 