Since $j$ increases by $1$ from term to term in the series, trivially, its increment is also $\Delta j= \left(j+1\right) - j = 1$, and it is correct to write
$$
\begin{aligned}
Z_\mathrm{tot} & = \sum_{j=0}^\infty(2j+1)e^{-j\left(j+1\right)\epsilon/kT} = \sum_{j=0}^\infty(2j+1)e^{-j\left(j+1\right)\epsilon/kT}\Delta j
\text{ .}
\end{aligned}
$$
Notice how this sum corresponds to a "left" rectangle integration rule of the function
$$
f(j)=(2j + 1)e^{-j\left(j+1\right)\epsilon/kT}\text{ ,}
$$
using subinvervals in $j$ of width $\Delta j$.
(Look at the yellow figure in the appended link, where the function value at the left limit of each subinterval is used as approximation for the function in each whole respective subinterval, as if the function was in fact rectangular.)
In the limit $\Delta j\to 0$, the sum would exactly match the respective integral, that is,
$$
\lim_{\Delta j\to 0} \sum_{j=0}^\infty(2j+1)e^{-j\left(j+1\right)\epsilon/kT} \Delta j=
\int_{j=0}^\infty(2j+1)e^{-j\left(j+1\right)\epsilon/kT}\,dj\text{ .}
$$
However, that is not the case here since $\Delta j=1$. Still, the sum may be a valid approximation to the integral if function $f$ does not vary too much in each interval $[j,\, j + \Delta j]$ where the function magnitude is significant (since where the function magnitude is not, that is, where $f(j)\approx 0$ regardless of the variation, the contribution of the respective terms in the sum would be close to null just like in the integral). In a circumstance of significant magnitude, one would need to have
$$
f(l)\approx f(j),\quad \text{for } l\in[j,\,j+\Delta j]\text{ .}
$$
Using a first order Tayler expansion around $l = j$ as an approximation to $f(l)$, one gets
$$
\begin{aligned}
f(l) & \approx f(j) + f'(j)\cdot(l - j)
\\
& \approx f(j) + (l - j)\cdot \left(2 - \frac{\epsilon}{kT}\left(2j+1\right)^2\right)e^{-j\left(j+1\right)\epsilon/kT}\text{ .}
\end{aligned}
$$
Thus,
$$
\begin{aligned}
\frac{\left|f(l)-f(j)\right|}{f(j)} &\le \left|\frac{2}{2j + 1} - \frac{\epsilon}{kT}\left(2j+1\right)\right|\underbrace{\Delta j}_{=1}
\\
& \le \left|\frac{2}{2j + 1} - \frac{\epsilon}{kT}\left(2j+1\right)\right|
\end{aligned}
$$
From this, one gets $\frac{\left|f(l)-f(j)\right|}{f(j)} \approx 0$ if $j\ll kT/\epsilon$, but $\frac{\left|f(l)-f(j)\right|}{f(j)} \gtrsim 1$ if $j\gtrsim kT/\epsilon$, meaning that the varition of function $f$ when compared to the function magnitude starts to be significant at $j\gtrsim kT/\epsilon$. Nevertheless, if the magnitude is negligible for such region, the approximation sum-to-integral would still be valid.
Well,
$$
f(j\sim \frac{kT}{\epsilon})\sim\frac{kT}{\epsilon}e^{-kT/\epsilon}\text{ .}
$$
Since $e^{-kT/\epsilon}$ grows much faster to $0$ than $\frac{kT}{\epsilon}$ to $\infty$, one has $f(j\sim \frac{kT}{\epsilon})\approx 0$. The function magnitude at $j\gtrsim \frac{kT}{\epsilon}$ would be much smaller than for previous values. Indeed, the maximum magnitude is mugh grater than 1:
$$
f_{\mathrm{max}}=f\left(\frac{1}{2}\sqrt{2\frac{kT}{\epsilon}}-\frac{1}{2}\right)=\underbrace{\sqrt{2 \frac{kT}{\epsilon}}}_{\gg 1}\cdot \underbrace{e^{\frac{1}{4}\frac{\epsilon}{kT} - \frac{1}{2}}}_{\sim 1}\gg 1\text{ .}
$$
One then concludes that under the condition $kT\gg\epsilon$, the approximation sum-to-integral is valid.