Burchnall-Chaundy polynomials and the Laurent phenomenon

The Burchnall-Chaundy polynomials $P_n(z)$ are determined by the differential recurrence relation $$P_{n+1}'(z)P_{n-1}(z)-P_{n+1}(z)P_{n-1}'(z)=P_n(z)^2$$ with $P_{-1}=P_0(z)=1.$ The fact that this recurrence relation has all solutions polynomial is not obvious and is similar to the integrality of Somos sequences and the Laurent phenomenon. We discuss this parallel in more detail and extend it to two difference equations $$Q_{n+1}(z+1)Q_{n-1}(z)-Q_{n+1}(z)Q_{n-1}(z+1)=Q_n(z)Q_n(z+1)$$ and $$R_{n+1}(z+1)R_{n-1}(z-1)-R_{n+1}(z-1)R_{n-1}(z+1)=R^2_n(z)$$ related to two different KdV-type reductions of the Hirota-Miwa and Dodgson octahedral equations. As a corollary we have a new form of the Burchnall-Chaundy polynomials in terms of the initial data $P_n(0)$, which is shown to be Laurent.