Actually $s$ doesn't need to systematically underestimate $\sigma$; this could happen even if that weren't true.
As it is, $s$ is biased for $\sigma$ (the fact that $s^2$ is unbiased for $\sigma^2$ means that $s$ will be biased for $\sigma$, due to Jensen's inequality*, but that's not the central thing going on there.
* Jensen's inequality
If $g$ is a convex function, $g\left(\text{E}[X]\right) \leq \text{E}\left[g(X)\right]$
with equality only if $X$ is constant or $g$ is linear.
Now $g(X)=-\sqrt{X}$ is convex,
so $-\sqrt{\text{E}[X]} < \text{E}(-\sqrt{X})$,
i.e.
$\sqrt{\text{E}[X]} > \text{E}(\sqrt{X})\,$, implying $\sigma>E(s)$ if the random variable $s$ is not a fixed constant.
Edit: a simpler demonstration not invoking Jensen --
Assume that the distribution of the underlying variable has $\sigma>0$.
Note that $\text{Var}(s) = E(s^2)-E(s)^2$ this variance will always be positive for $\sigma>0$.
Hence $E(s)^2 = E(s^2)-\text{Var}(s) < \sigma^2$, so $E(s)<\sigma$.
So what is the main issue?
Let $Z=\frac{\overline{X} - \mu}{\frac{\sigma}{\sqrt{n}}}$
Note that you're dealing with $t=Z\cdot\frac{\sigma}{s}$.
That inversion of $s$ is important. So the effect on the variance it's not whether $s$ is smaller than $\sigma$ on average (though it is, very slightly), but whether $1/s$ is larger than $1/\sigma$ on average (and those two things are NOT the same thing).
And it is larger, to a greater extent than its inverse is smaller.
Which is to say $E(1/X)\neq 1/E(X)$; in fact, from Jensen's inequality:
$g(X) = 1/x$ is convex, so if $X$ is not constant,
$1/\left(\text{E}[X]\right) < \text{E}\left[1/X\right]$
So consider, for example, normal samples of size 10; $s$ is about 2.7% smaller than $\sigma$ on average, but $1/s$ is about 9.4% larger than $1/\sigma$ on average. So even if at n=10 we made our estimate of $\sigma$ 2.7-something percent larger** so that $E(\widehat\sigma)=\sigma$, the corresponding $t=Z\cdot\frac{\sigma}{\widehat\sigma}$ would not have unit variance - it would still be a fair bit larger than 1.
**(at other $n$ the adjustment would be different of course)
Since the t-distribution is like the standard normal distribution but with a higher variance (smaller peak and fatter tails)
If you adjust for the difference in spread, the peak is higher.
Why does the t-distribution become more normal as sample size increases?
The standard normal distribution vs the t-distribution
One problem I see is that in Gamma regression the shape parameter is assumed constant, and the scale parameter is the one that changes.
However, if your underlying variables are $N(\mu_i,\sigma^2_i)$ from which you draw a sample of size $n_i$ and compute the sample variance $s^2_i$, then $(n_i-1)s^2/\sigma^2_i$ will be $\chi^2_{n-1}$, and a $\chi^2_{\nu}$ is a $\text{Gamma}(\nu/2,2)$ (shape-scale parameterization). So for your data it seems that the shape parameter varies with $n-1$, but not the scale parameter.
Note that as $n-1$ increases the shape of the chi-square becomes less skew, but the shape in the Gamma GLM modelled as you suggest does not become less skew as $n-1$ increases.
You might even be better working with an overdispersed Poisson model; at least it will become less skew as $n-1$ increases, unlike the Gamma model.
Even better perhaps would be to model the likelihood more directly.
Best Answer
Poisson regression finds a value $\hat{\beta}$ maximizing the likelihood of the data. For any value of $x$, you would then suppose $Y$ has a Poisson($\exp(x \hat{\beta})$) distribution. The standard deviation of that distribution equals $\exp(x \hat{\beta}/2)$. This appears to be what you mean by $\sqrt{\widehat{\mu}}$.
There are, of course, other ways to estimate the standard deviation of $Y|x$. However, staying within the context of Poisson regression, $\exp(x \hat{\beta}/2)$ is the ML estimator of SD($Y|x$) for the simple reason that the ML estimator of a function of the parameters is the same function of the ML estimator of those parameters. The function in this case is the one sending $\hat{\beta}$ to $\exp(x \hat{\beta}/2)$ (for any fixed value of $x$). This theorem will appear in any full account of maximum likelihood estimation. Its proof is straightforward. Conceptually, the function is a way to re-express the parameters, but re-expressing them doesn't change the fact that they maximize (or fail to maximize, depending on their values) the likelihood.