Computational studies of transcription factor binding to DNA are usually based on a simple matrix model (PWM) of sequence-dependent binding energy. For some transcription factors, the best matrix we can construct is alarmingly non-specific, predicting many binding sites that are not known to be functional. If these sites are in fact spurious, the validity of the matrix model as a physical description of binding specificity is called into question. This is unfortunate, since the model is a very convenient framework for physical thinking about gene regulatory networks. This has motivated us to develop a species comparison approach to assessing the functionality of populations of such sites: we compare putative sites with aligned sequences in other genomes and ask whether the probability of mutation depends on position within the site. We find a very distinct `footprint' in the mutation probability which has the feature that mutation is most strongly suppressed where it would have the biggest effect on the binding energy as assessed by the matrix model. This suggests that many of the apparently spurious sites are functional and also that something related to the matrix model estimate of site binding energy is what is conserved between species. Apart from theoretical reassurance, this analysis casts up new binding sites and new regulated genes which may be of interest to experiment.