Let for α, β, γ > −1, α+γ+ 3 2 > 0, β+γ+ 3 2 > 0 and n ⩾ k ⩾ 0 the orthogonal polynomials p n, k α, β, γ ( u, v) be defined as polynomials in u and v with “highest” term u n− k v k which are obtained by orthogonalization of the sequence 1, u, v, u 2, uv, v 2, u 3, u 2 v, … with respect to the weight function (1− u+ v) α (1+ u+ v) β ( u 2−4 v) γ on the region bounded by the lines I− u+ v=0 and 1+ u+ v−0 and by the parabola u 2−4 v=0. Two explicit linear partial differential operators D 1 α, β, γ and D 2 α, β, γ of orders two and four, respectively, are obtained such that the polynomials p n, k α, β, γ ( u, v) are eigenfunctions of D 1 α, β, γ and D 2 α, β, γ . It is proved that if a differential operator D has the polynomials p n, k α, β, γ ( u, v) as eigenfunctions then D can be expressed in one and only one way as a polynomial in D 1 α, β, γ and D 2 α, β, γ . The special case γ=− 1 2 can be reduced to Jacobi polynomials by the identity p n, k α,β,− 1 2 (x+y, xy)= const. (P n (α,β)(x)P k (α,β)(y)+P k (α,β)(x)P n (α,β)(y)) . For certain values of α, β, γ and in terms of the coordinates s, t, where u=cos s+cos t, v=cos s cos t, the operator D 1 α, β, γ is the radial part of the Laplace-Beltrami operator on certain compact Riemannian symmetric spaces of rank two.