Summary: The integral $\int_0^\infty t^{\frac{s-1}{2}} \operatorname{Tr}[ D^u (\exp(-t D^2) ] \, dt$ has a simple pole at $s=0$ with residue $2 C_{-1/2}$, from which your expression for $\eta(0)$ follows. To compute the residue, we compute the integral explicitly (up to a holomorphic remainder) near $t=0$ for $\operatorname{Re} s$ large, then analytically continue. Each term in the $t \to 0^+$ asymptotic expansion of the heat trace corresponds to a pole in $s$ of the integral.
To explain, I'm going to phrase this as a more general property of zeta functions defined by a heat trace via the Mellin transform; then at the end we'll return to the eta function. I think these ideas are first due to Seeley in the paper Complex powers of an elliptic operator. A good reference is Gilkey's book Invariance theory, the heat equation, and the Atiyah-Singer index theorem. But I will try to write out all the details here.
Let $h(t)$ ($t > 0$) be a smooth function that satisfies these two conditions:
(A): $h(t)$ is $O(\exp(-\epsilon t ))$ as $t \to \infty$ (for some $\epsilon > 0$).
(B): $h(t) = \displaystyle \sum_{k=0}^K C_{p_k} t^{p_k} + R_K(t)$, where $R_K(t)$ is $O(t^{p_{K+1}})$ as $t \to 0^+$ (for some increasing sequence of powers $p_0 < p_1 < p_2 < \dots \to \infty$).
Example/remark: The original question is interested in the case when $h(t) = \operatorname{Tr}[ D^u \exp(-t D^2) ]$, which satisfies those conditions with $p_k = -n/2 + k$, as you said. Note that $\operatorname{Tr}[ P \exp(-t L) ]$ also satisfies the conditions for more general operators $P$ and $L$, where $L$ is a strictly positive elliptic differential operator. (Or if $L$ is nonnegative but has a nontrivial kernel, we can project away from the kernel to get the exponential decay as $t\to \infty$.)
The zeta function associated to $h(t)$ is defined using the Mellin transform. The integral in the following definition converges for $s \in \mathbb{C}$ with $\operatorname{Re} s > -p_0$ (see below):
$$ \zeta(s) := \frac{1}{\Gamma(s)} \int_0^\infty t^{s-1} h(t) \, dt.$$
Note that $\frac{1}{\Gamma(s)}$ is entire (i.e., holomorphic everywhere) since the gamma function $\Gamma(s)$ is nonvanishing. Thus to study poles of $\zeta(s)$, we really just need to study the integral
$$ I(s) := \int_0^\infty t^{s-1} h(t) \, dt = I_A(s) + I_B(s),$$
where the two pieces are defined by
$$ I_A(s) := \int_1^\infty t^{s-1} h(t) \, dt, ~~~~ I_B(s) := \int_0^1 t^{s-1} h(t) \, dt.$$
We'll deal with the two pieces using the respective conditions above:
(A): Since $h(t)$ is $O(\exp(-\epsilon t ))$ as $t \to \infty$, the integral defining $I_A(s)$ converges for all $s \in \mathbb{C}$, and $I_A(s)$ is entire (so we can ignore it for the sake of studying poles).
(B): Using the asymptotic expansion, for each $K$ (which defines the point at which we truncate the asymptotic expansion), we have
\begin{align*}
I_B(s) &= \sum_{k=0}^K C_{p_k} \int_0^1 t^{s-1 + p_k} \, dt + \int_0^1 t^{s-1} R_K(t) \, dt \\
&= \sum_{k=0}^K C_{p_k} \frac{1}{s+p_k} + \text{(holomorphic for $\operatorname{Re} s > -p_{K+1}$ by the estimate on $R_K(t)$)}
\end{align*}
This shows that $I_B(s)$ is holomorphic for $\operatorname{Re} s > -p_0$
(i.e., $\operatorname{Re} s > \frac{n}{2}$ in your case), and furthermore
that $I_B(s)$ possesses an analytic continuation to the half-plane where
$\operatorname{Re} s > -p_{K+1}$ with (for $k=0, 1, 2, \dots, K$) a simple pole at $s = -p_k$, where the residue is $C_{p_k}$.
Putting all that together, we have that $\zeta(s) = \frac{1}{\Gamma(s)}(I_A(s) + I_B(s))$ has (for $k=0, 1, 2, \dots$) a simple pole at $s= -p_k$, where the residue is $\frac{C_{p_k}}{\Gamma(-p_k)}$. In particular, we will need that the residue of $\zeta(s)$ at $s = \frac{1}{2}$ is $\frac{C_{-1/2}}{\Gamma(1/2)} = \frac{C_{-1/2}}{\sqrt{\pi}}$. (I looked up $\Gamma(1/2)$ on Wikipedia.)
Now finally we can return to your question about the eta function. If I understand correctly, your definition of $\eta(z)$ is (this appears to be a derivative of some other eta function)
$$ \eta(z) = -z \zeta \left( \frac{z+1}{2} \right).$$
(I hope the new variable $z$ helps avoid confusion.) Since $\zeta \left( \frac{z+1}{2} \right)$ has only simple poles, the factor of $-z$ gives that $\eta(0)$ is $-1$ times the residue of $\zeta \left( \frac{z+1}{2} \right)$ at $z=0$, which is $2$ times the residue of $\zeta(s)$ at $s = \frac{1}{2}$, which we computed above to be $\frac{C_{-1/2}}{\sqrt{\pi}}$. Thus
$$\eta(0) = \frac{-2 C_{-1/2}}{\sqrt{\pi}}.$$
