Your (unbounded) operator is in the limit circle case at $+\infty$. This means that self-adjoint realizations are obtained only when a boundary condition at $\infty$ is imposed, and the eigenvalues will depend on this boundary condition. (In fact, for any $\lambda\in\mathbb R$, there is exactly one boundary condition that makes this $\lambda$ an eigenvalue.)
To see this, let's make the (obvious) try $u_n=q^{-2n}g_n$. Then $Tg=\lambda g$ (viewed as a difference equation) is equivalent to
$$
q^{1/2}u_{n+1}+q^{-3/2}u_{n-1}+\frac{1+q}{q}u_n = \lambda q^{2n} u_n .
$$
The coefficient $\lambda q^{2n}$ is summable on $n\ge 1$, so can be ignored as far as the asymptotic behavior of $u_n$ is concerned, and we are in control. (To prove this carefully, we would have to refer to a discrete version of Levinson's theorem, but in fact this is unnecessary for our current purposes since it's enough to observe that all solutions are square summable when $\lambda =0$.)
Two basis solutions for $\lambda =0$ are $u_n=(-q)^{-n/2}$ and $u_n=(-q)^{-3n/2}$. Going back to $g$, we thus have the two solutions $|g_n|=q^{3n/2}$ and $|g_n|= q^{n/2}$, both of which are exponentially decaying and in $\ell^2(\mathbb Z_+)$. All solutions are square integrable, as claimed.
We can also observe that the coefficients of the left half line operator $T_-$ go to zero, so $\sigma_{ess}(T_-)=\{ 0\}$. This shows that (for any boundary condition at $\infty$) $\sigma_{ess}(T)=\{ 0\}$, so the eigenvalues accumulate at $0$ and the spectrum is purely discrete everywhere else.