I think there are several things you could do and it is hard to point out all of them here.
In the first place I would look at the phylogenetic relationship of the strains, meaning 16S, core genome phylogenesis, SNP.
This imply that you would need to calculate core and pan genome as well. This analysis will also give you good information on the gene uniquely present in specific strains. Looking into them can maybe give you some insight into potential uniq features, e.g. virulence.
Then I would look at those system involved in virulence: secretion systems (II, III, IV, VI, VII) and also their effectors. This are crucial in establishing a relationship with other bacteria and eukaryotes. It could be that the virulent strains have specific effectors that the other do not have.
There are then specific db for virulence factors, you could take this proteins and blast them against your genomes. I would not recommend this approach thought, because you get many false hits.
Finally, you could see the dNdS among the strains, maybe the virulence strains have specific genes under positive selection compared to the other. You can apply this only to orthologous genes and you should be careful if the strain are to close to each other. Often strains that are too closely related show higher dNdS than aspected. It seems that this is not a real effect of positive selection, but rather a bias due to the little divergence time between the strains.
Hope this helps
I have toggled the two most upvoted answers as the accepted answers as OP has not been seen in 5+ years.
Multiple genome alignments (e.g. with mauve) may give you insights about syntenic regions and recombination events which are mostly related the pathogenic potential.