Inspired by your comparison with the Hurwitz zeta function, which has the integral representation $$ \zeta (s,a)={\frac {1}{\Gamma (s)}}\int _{0}^{\infty }{\frac {x^{s-1}e^{-ax}}{1-e^{-x}}}\mathrm{d}x \:, $$ I tried to come up with an integral representation for your function $f$. The factor $s$ that comes out when taking the $q$ derivative of $\zeta(q,s)$$\zeta(s,q)$ arises from the $\Gamma(s)$ in the denominator. To get your factor $\sqrt{s(2+s)}$, we need to change this expression into something of the form $$ f(q,s) = \frac{1}{\sqrt{\Gamma(s+1)\Gamma(s+3)}} \int_0^\infty x^{s-1} e^{-q x} G(x) \mathrm{d} x \:, $$ with an unknown function $G$. This expression satisfies your PDE. The function $G$ can be determined from your condition $f(0,s) = 1$, which actually gives the Mellin transform of $G$, $$ \int_0^\infty x^{s-1} G(x) \mathrm{d}x = \sqrt{\Gamma(s+1)\Gamma(s+3)} = \Gamma(s) s\sqrt{(s+2)(s+1)} \:. $$ In this form, you can obtain $G$ as a power series via Ramanujan's master theorem, which gives $$ G(x) = -\sum_{k=3}^\infty \frac{(-x)^k}{(k-1)!}\sqrt{(k-1)(k-2)} \:. $$
