It looks like the tentative conjecture is true -- $\zeta_\theta(s)$ extends to a meromorphic function on $\mathbb{C}$ for all real algebraic numbers. In fact, it will extend to a meromorphic function whenever $\theta$ cannot be approximated by rationals (in particular, if $\theta$ has irrationality measure 2 which, by the Thue-Siegel-Roth theorem, includes all algebraic irrationals). Conversely, if an irrational $\theta$ can be approximated too well by rationals then there will not be a meromorphic extension, which will be the case outside of a meagre subset of $\mathbb{R}^+$. *[I think the following is all correct, but I would be more confident if someone double checks it...]*

First, let $\lbrace x\rbrace = x - \lfloor x\rfloor$ denote the fractional part of $x$, and write
$$
\tilde\zeta_\theta(s)=\sum_{n\ge1}\frac{\lbrace n\theta\rbrace-1/2}{n^s},
$$
which converges on $\Re(s) > 1$ (at the very least).

Then, the function $\zeta_\theta$ defined in the question is
$$
\zeta_\theta(s)=\theta\zeta(s-1)-\frac12\zeta(s)-\tilde\zeta_\theta(s),
$$
where $\zeta$ is the Riemann zeta function. So, at least, we have an extension to $\Re(s) > 1$. In the case where $\theta$ is rational then $\lbrace n\theta\rbrace-1/2$ is periodic in $n$, so $\tilde\zeta_\theta$ is a linear combination of Dirichlet L-functions and has a meromorphic extension to $\mathbb{C}$. Henceforth, I will only consider irrational $\theta$. Setting $F(x)=\sum_{1\le n\le x}(\lbrace n\theta\rbrace-1/2)$, the equidistribution theorem states that $F(x)/x\to0$ as $x\to\infty$. Rearranginging the expression for $\tilde\zeta_\theta$ gives
$$
\tilde\zeta_\theta(s)=s\sum_{n\ge1}\frac{F(n)}{n}n^{-s}\left(\frac ns(1-(1+1/n)^{-s})\right).
$$
If $\lvert F(n)/n\rvert$ is bounded by some $\epsilon > 0$ then this expression will be bounded by $\epsilon\zeta(s)$ for all $\Re(s) > 1$. So, by equidistribution, it follows that $(s-1)\tilde\zeta_\theta(s)\to0$ as $s\to1$ (on $s > 1$). If a meromorphic extension existed in a neighbourhood of 1 then this would imply that $\tilde\zeta_\theta$ is bounded near 1. However, if $\theta=p/q$ is rational ($p,q$ coprime) then $F(x)/x\to-1/(2q)$ so, if $x$ is too closely approximated by rationals then there will be a slight bias in $F(x)/x$ and $\tilde\zeta_\theta(s)$ will not be bounded close to 1. I think that the existence of rational approximations $p_n/q_n$ with $q_n^2\log\lvert\theta-p_n/q_n\rvert\to-\infty$ is enough for this to happen.

Now, suppose that $\theta$ has irrationality measure 2, so that there are only finitely many rational approximations $\lvert\theta-p/q\rvert\le q^{-2-\epsilon}$ for each $\epsilon > 0$. Equivalently, there are constants $C_\epsilon > 0$ such that $n^{1+\epsilon}\lvert n\theta-m\rvert\ge C_\epsilon$ for all positive integer $n$ and integer $m$. Then, for each $\epsilon > 0$, the discrepancy of the set $\left\lbrace \lbrace k\theta\rbrace\colon 1\le k\le n\right\rbrace$ is bounded by $Kn^{-1+\epsilon}$ (for a constant $K$, depending on $\epsilon > 0$). This means that $F(n)/n$ is bounded by $2Kn^{-1+\epsilon}$ and, hence, the expression above for $\tilde\zeta(s)$ converges for all $\Re(s) > \epsilon$. So, $\tilde\zeta(s)$ is an analytic function on $\Re(s) > 0$. Therefore, $\zeta_\theta(s)$ is meromorphic in this region with simple poles at $s=2$ and $s=1$.

We can now extend to the whole of $\mathbb{C}$ via a functional equation. Hurwitz's formula gives, for $\Re(s) > 1$,
$$
\begin{array}{l}
&\sum_{n\ge1}e^{2\pi inx}n^{-s}=(2\pi)^{s-1}\Gamma(1-s)\left(e^{i\pi(1-s)/2}\zeta(1-s,x)+e^{i\pi(s-1)/2}\zeta(1-s,1-x)\right)
\end{array}
$$
where $\zeta(s,x)=\sum_{n\ge1}(n+x)^{-s}$ is the Hurwitz Zeta function.
Plugging in the Fourier series $\lbrace x\rbrace=1/2+\sum_{k\not=0}(i/(2\pi k))e^{2\pi ikx}$, applying Hurwitz's formula, and being careless about when the sums converge/commute (it can be made more rigorous) gives,
$$
\begin{array}{rl}
\tilde\zeta_\theta(s)&=\sum_{n\ge1}\sum_{k\not=0}\frac{ie^{2\pi ikn\theta}}{2\pi kn^s}\cr
&=\sum_{k\not=0}\frac{i}{2\pi k}\sum_{n\ge1}\frac{e^{2\pi in\lbrace k\theta\rbrace}}{n^s}\cr
&=(2\pi)^{s-2}\Gamma(1-s)\sum_{k\not=0}\frac{i}{k}\left(e^{i\pi(1-s)/2}\zeta(1-s,\lbrace k\theta\rbrace)+e^{i\pi(s-1)/2}\zeta(1-s,1-\lbrace k\theta\rbrace)\right)\cr
&=(2\pi)^{s-2}\Gamma(1-s)\sum_{k\not=0}\frac{i}{k}\left(e^{i\pi(1-s)/2}\zeta(1-s,\lbrace k\theta\rbrace)-e^{i\pi(s-1)/2}\zeta(1-s,\lbrace k\theta\rbrace)\right)\cr
&=2\sin((s-1)\pi/2)(2\pi)^{s-2}\Gamma(1-s)\sum_{k\not=0}k^{-1}\zeta(1-s,\lbrace k\theta\rbrace)\cr
&=2\sin((s-1)\pi/2)(2\pi)^{s-2}\Gamma(1-s)\sum_{n\ge1}\sum_{k\not=0}k^{-1}(n+\lbrace k\theta\rbrace)^{s-1}.
\end{array}
$$
Setting
$$
\begin{array}{l}
&g_{n,x,s}=(1+\lbrace x\rbrace/n)^{-s}-(1+(1-\lbrace x\rbrace)/n)^{-s},\cr &G_{n,x,s}=\sum_{1\le k\le x}g_{n,k\theta,s},
\end{array}
$$
we have
$$
\begin{array}{rl}
\tilde\zeta_\theta(s)&=2\sin((s-1)\pi/2)(2\pi)^{s-2}\Gamma(1-s)\sum_{n\ge1}\sum_{k\ge1}k^{-1}g_{n,k\theta,1-s}n^{s-1}\cr
&=2\sin((s-1)\pi/2)(2\pi)^{s-2}\Gamma(1-s)\sum_{n\ge1}\sum_{k\ge1}\frac{G_{n,k,1-s}}{k(k+1)}n^{s-1}.
\end{array}
$$
Now, as the variation of $g_{n,x,s}$ on $0\le x\le1$ is bounded by $\lvert s\rvert/n$ for $\Re(s) > 0$, $G_{n,k,s}$ is of size $O(\lvert s\rvert k^\epsilon /n)$ and, hence, the sum $\sum_{k\ge1}G_{n,k,s}/(k(k+1))$ is of size $O(\lvert s\rvert/n)$. Therefore, the final expression above for $\tilde\zeta_\theta(s)$ is a sum over $n\ge1$ of a term of size $O(\lvert 1-s\rvert n^{s-2})$, which converges (uniformly on compacts) to an analytic function on $\Re(s) < 1$, giving the desired extension of $\tilde\zeta_\theta$ to $\mathbb{C}$. Note that I have arrived at $\tilde\zeta_\theta$ entire, so that $\zeta_\theta$ only has poles at $s=1$ and $s=2$. There are no poles on the imaginary axis as suggested in the question.