The intuition about the "plus" signs related to the variance (from the fact that even when we calculate the variance of a difference of independent random variables, we add their variances) is correct but fatally incomplete: if the random variables involved are not independent, then covariances are also involved -and covariances may be negative. There exists an expression that is almost like the expression in the question was thought that it "should" be by the OP (and me), and it is the variance of the prediction error, denote it $e^0 = y^0 - \hat y^0$, where $y^0 = \beta_0+\beta_1x^0+u^0$:
$$\text{Var}(e^0) = \sigma^2\cdot \left(1 + \frac 1n + \frac {(x^0-\bar x)^2}{S_{xx}}\right)$$
The critical difference between the variance of the prediction error and the variance of the estimation error (i.e. of the residual), is that the error term of the predicted observation is not correlated with the estimator, since the value $y^0$ was not used in constructing the estimator and calculating the estimates, being an out-of-sample value.
The algebra for both proceeds in exactly the same way up to a point (using $^0$ instead of $_i$), but then diverges. Specifically:
In the simple linear regression $y_i = \beta_0 + \beta_1x_i + u_i$, $\text{Var}(u_i)=\sigma^2$, the variance of the estimator $\hat \beta = (\hat \beta_0, \hat \beta_1)'$ is still
$$\text{Var}(\hat \beta) = \sigma^2 \left(\mathbf X' \mathbf X\right)^{-1}$$
We have
$$\mathbf X' \mathbf X= \left[ \begin{matrix}
n & \sum x_i\\
\sum x_i & \sum x_i^2 \end{matrix}\right]$$
and so
$$\left(\mathbf X' \mathbf X\right)^{-1}= \left[ \begin{matrix}
\sum x_i^2 & -\sum x_i\\
-\sum x_i & n \end{matrix}\right]\cdot \left[n\sum x_i^2-\left(\sum x_i\right)^2\right]^{-1}$$
We have
$$\left[n\sum x_i^2-\left(\sum x_i\right)^2\right] = \left[n\sum x_i^2-n^2\bar x^2\right] = n\left[\sum x_i^2-n\bar x^2\right] \\= n\sum (x_i^2-\bar x^2) \equiv nS_{xx}$$
So
$$\left(\mathbf X' \mathbf X\right)^{-1}= \left[ \begin{matrix}
(1/n)\sum x_i^2 & -\bar x\\
-\bar x & 1 \end{matrix}\right]\cdot (1/S_{xx})$$
which means that
$$\text{Var}(\hat \beta_0) = \sigma^2\left(\frac 1n\sum x_i^2\right)\cdot \ (1/S_{xx}) = \frac {\sigma^2}{n}\frac{S_{xx}+n\bar x^2} {S_{xx}} = \sigma^2\left(\frac 1n + \frac{\bar x^2} {S_{xx}}\right) $$
$$\text{Var}(\hat \beta_1) = \sigma^2(1/S_{xx}) $$
$$\text{Cov}(\hat \beta_0,\hat \beta_1) = -\sigma^2(\bar x/S_{xx}) $$
The $i$-th residual is defined as
$$\hat u_i = y_i - \hat y_i = (\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i +u_i$$
The actual coefficients are treated as constants, the regressor is fixed (or conditional on it), and has zero covariance with the error term, but the estimators are correlated with the error term, because the estimators contain the dependent variable, and the dependent variable contains the error term. So we have
$$\text{Var}(\hat u_i) = \Big[\text{Var}(u_i)+\text{Var}(\hat \beta_0)+x_i^2\text{Var}(\hat \beta_1)+2x_i\text{Cov}(\hat \beta_0,\hat \beta_1)\Big] + 2\text{Cov}([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i],u_i) $$
$$=\Big[\sigma^2 + \sigma^2\left(\frac 1n + \frac{\bar x^2} {S_{xx}}\right) + x_i^2\sigma^2(1/S_{xx}) +2\text{Cov}([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i],u_i)$$
Pack it up a bit to obtain
$$\text{Var}(\hat u_i)=\left[\sigma^2\cdot \left(1 + \frac 1n + \frac {(x_i-\bar x)^2}{S_{xx}}\right)\right]+ 2\text{Cov}([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i],u_i)$$
The term in the big parenthesis has exactly the same structure with the variance of the prediction error, with the only change being that instead of $x_i$ we will have $x^0$ (and the variance will be that of $e^0$ and not of $\hat u_i$). The last covariance term is zero for the prediction error because $y^0$ and hence $u^0$ is not included in the estimators, but not zero for the estimation error because $y_i$ and hence $u_i$ is part of the sample and so it is included in the estimator. We have
$$2\text{Cov}([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i],u_i) = 2E\left([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i]u_i\right)$$
$$=-2E\left(\hat \beta_0u_i\right)-2x_iE\left(\hat \beta_1u_i\right) = -2E\left([\bar y -\hat \beta_1 \bar x]u_i\right)-2x_iE\left(\hat \beta_1u_i\right)$$
the last substitution from how $\hat \beta_0$ is calculated. Continuing,
$$...=-2E(\bar yu_i) -2(x_i-\bar x)E\left(\hat \beta_1u_i\right) = -2\frac {\sigma^2}{n} -2(x_i-\bar x)E\left[\frac {\sum(x_i-\bar x)(y_i-\bar y)}{S_{xx}}u_i\right]$$
$$=-2\frac {\sigma^2}{n} -2\frac {(x_i-\bar x)}{S_{xx}}\left[ \sum(x_i-\bar x)E(y_iu_i-\bar yu_i)\right]$$
$$=-2\frac {\sigma^2}{n} -2\frac {(x_i-\bar x)}{S_{xx}}\left[ -\frac {\sigma^2}{n}\sum_{j\neq i}(x_j-\bar x) + (x_i-\bar x)\sigma^2(1-\frac 1n)\right]$$
$$=-2\frac {\sigma^2}{n}-2\frac {(x_i-\bar x)}{S_{xx}}\left[ -\frac {\sigma^2}{n}\sum(x_i-\bar x) + (x_i-\bar x)\sigma^2\right]$$
$$=-2\frac {\sigma^2}{n}-2\frac {(x_i-\bar x)}{S_{xx}}\left[ 0 + (x_i-\bar x)\sigma^2\right] = -2\frac {\sigma^2}{n}-2\sigma^2\frac {(x_i-\bar x)^2}{S_{xx}}$$
Inserting this into the expression for the variance of the residual, we obtain
$$\text{Var}(\hat u_i)=\sigma^2\cdot \left(1 - \frac 1n - \frac {(x_i-\bar x)^2}{S_{xx}}\right)$$
So hats off to the text the OP is using.
(I have skipped some algebraic manipulations, no wonder OLS algebra is taught less and less these days...)
SOME INTUITION
So it appears that what works "against" us (larger variance) when predicting, works "for us" (lower variance) when estimating. This is a good starting point for one to ponder why an excellent fit may be a bad sign for the prediction abilities of the model (however counter-intuitive this may sound...).
The fact that we are estimating the expected value of the regressor, decreases the variance by $1/n$. Why? because by estimating, we "close our eyes" to some error-variability existing in the sample,since we essentially estimating an expected value. Moreover, the larger the deviation of an observation of a regressor from the regressor's sample mean, the smaller the variance of the residual associated with this observation will be... the more deviant the observation, the less deviant its residual... It is variability of the regressors that works for us, by "taking the place" of the unknown error-variability.
But that's good for estimation. For prediction, the same things turn against us: now, by not taking into account, however imperfectly, the variability in $y^0$ (since we want to predict it), our imperfect estimators obtained from the sample show their weaknesses: we estimated the sample mean, we don't know the true expected value -the variance increases. We have an $x^0$ that is far away from the sample mean as calculated from the other observations -too bad, our prediction error variance gets another boost, because the predicted $\hat y^0$ will tend to go astray... in more scientific language "optimal predictors in the sense of reduced prediction error variance, represent a shrinkage towards the mean of the variable under prediction". We do not try to replicate the dependent variable's variability -we just try to stay "close to the average".