By analogy with the elementary example of the harmonic oscillator, of Calogero-Sutherland-Moser integrable models or the open Toda chain for which a (q-)hypergeometric formulation of the Hamiltonian's eigenfunctions i known (in terms of one-variable Hermite polynomials, multivariable Macdonald-Koornwinder or limiting cases (q-)Wittaker or Hall-Littlewood), I will show that the bispectral family of multivariable Gasper-Rahman polynomials (generalizing those of Askey-Wilson) provide a q-hypergeometric basis for the Hamiltonian's eigenfunctions of the class of quantum integrable models generated from the q-Onsager algebra (ex: Ising, superintegrable chiral Potts, XY, open XXZ,...). In a first part, the structure of Gasper-Rahman polynomials and their properties (orthogonality, bispectrality) will be recalled. In a second part, it will be shown that such polynomials provide infinite or finite modules of the q-Onsager algebra. In a third part, the action of elements of the Abelian subalgebra (generating the integrals of motion) in the polynomial basis will be considered. In particular, the conditions under which invariant subspaces exist will be identified, giving a new interpretation of Nepomechie's constraint. Consequences and perspectives for above mentionned integrable models as well as the interest for other formulations (Bethe ansatz, SOV) will be briefly described. \ \ Work done in collaboration with X. Martin (Tours).