Should I include solvation in quantum chemical calculations of a simplified protein-ligand interactions?
This is difficult to answer without full knowledge of the system, but some tips:

1) The overall energy of the system is arguably irrelevant.  Rather, the relative difference between different conformations, and perhaps, different ligands is of utility.

2) The level of theory and basis set your using is quite low.  If you have the time and the resources, scale it up.  For biomolecules I've found the minnesota functionals to be most accurate (if insisting on DFT) but I still believe MP2 is king.

3) Referring back to #1, if you solvate, then your comparison should too be solvated.  Perhaps your approach would be to run a solvated and unsolvated model for both your ligand, and perhaps, another known ligand. This could too help convince you of the utility of the solvation model.