The Heine-Stieltjes theorem describes the polynomial solutions, (v,f) such that T(f)=vf, to specific second-order differential operators, T, with polynomial coefficients. We extend the theorem to concern all (nondegenerate) differential operators preserving the property of having only real zeros, thus solving a conjecture of B. Shapiro. The new methods developed are used to describe intricate interlacing relations between the zeros of different pairs of solutions. This extends recent results of Bourget, McMillen and Vargas for the Heun equation and answers their question of how to generalize their results to higher degrees. Many of the results are new even for the classical case.