The question whether or not to exponentiate, is a question of what scale is employed for the normal approximation for the confidence intervals.
In negative binomial regressions you usually use a log-link, i.e., $\log(\mu) = \eta = x^\top \beta$, where $\mu$ is the expectation of the response and $\eta = x^\top \beta$ is the corresponding linear predictor. Given fitted regression coefficients $\hat \beta$ one can easily obtain the fitted linear predictor $\hat \eta = x^\top \hat \beta$ (called type = "link" for glm.nb()) and the corresponding expected response $\hat \mu = \exp(\hat \eta)$ (called type = "response" for glm.nb()).
However, for obtaining the corresponding confidence intervals it makes a difference whether the normal approximation is applied for $\hat \eta$ or for $\hat \mu$. Asymptotically, both versions agree but in finite samples they are different. The predict() method for glm.nb() can return standard errors on either scale but some might argue that it is more natural to apply the normal approximation on the linear predictor scale. The reason is simply that the linear predictor lives on the real line (like the normal distribution) and subsequently applying the inverse link $\exp(\cdot)$ assures positivity of the interval. So I would guess that the blog post you are referring to took this approach.
In practice, however, the difference may be rather small. As a simple reproducible example (that is a really non-sensical but requires no special packages) consider:
library("MASS")
m <- glm.nb(dist ~ speed, data = cars)
So we model how stopping distance increases with speed. (It's not really count data...but it's easy to visualize and does not require any further R packages here.) Then we can add the fitted values and standard errors on the response ($\mu$) scale and on the link ($\eta$) scale:
cars <- cbind(cars, "resp" = predict(m, type = "response", se.fit = TRUE)[1:2])
cars <- cbind(cars, "link" = predict(m, type = "link", se.fit = TRUE)[1:2])
Plotting the data, the fitted expectation, and the confidence intervals from the response scale (in red) yields the figure on the left below.
This is obtained with the usual approximation adding 2 times the standard errors:
plot(dist ~ speed, data = cars)
lines(resp.fit ~ speed, data = cars, col = 2)
lines(resp.fit + 2 * resp.se.fit ~ speed, data = cars, col = 2)
lines(resp.fit - 2 * resp.se.fit ~ speed, data = cars, col = 2)
Adding also the exponentiated fitted values and confidence interval from the link scale (in blue) then yields the plot on the right.
lines(exp(link.fit) ~ speed, data = cars, col = 4)
lines(exp(link.fit + 2 * link.se.fit) ~ speed, data = cars, col = 4)
lines(exp(link.fit - 2 * link.se.fit) ~ speed, data = cars, col = 4)
You can see that the fitted responses exactly coincide but the confidence intervals are slightly different, especially for higher values. Of course, the differences may or may not be larger on other data - but generally they will decrease with the sample size.
