The reason is that the compounded return of an investment has a "drag" when the intermediate returns are more volatile. Check for instance the compounded return of two return series with the same arithmetic mean, but different variance. You will see that the series with a lower variance has a higher compounded return (and geometric mean). The reason is that if an investment looses 50% of its value, it has to make a return of +100% (very different from +50%) to come back to the initial value.
The value of a bank account with no risk paying a continuously compounded interest rate $r$ is $S_0e^{r t}$. However, the GBM $S(t)$ is a model for the value of a risky/stochastic investment.
As the exact solution for $S(t)$ shows, the correct interpretation of the value of $S(t)$ as paying a continuously compounded rate of return is with $\gamma = r - \frac{1}{2}\sigma^2$ (also called growth rate) and not $r$, because $$\underset{t\rightarrow \infty}{\lim} S(t)=S_0e^{\gamma t}$$ Although the mean of $S(t)$ is $E[S(t)]=S_0e^{r t}$, the median of its distribution is $S_0e^{\gamma t}$.
While the intuition of an expected value as a proxy for the long-term value is appropriate for random variables with broadly symmetric distributions, such as returns, $dS(t)/S(t)$ (and other stationary series), in the case of the value of $S(t)$ this intuition is flawed.
To see this, notice that if $\gamma>0$ then $\underset{t\rightarrow \infty}{\lim} S(t)=\infty$, and if $\gamma<0$ then $\underset{t\rightarrow \infty}{\lim} S(t)=0$. On the other hand, $r>0$ does not guarantee an increasing value of $S$ in the long run, because if $\sigma$ is large enough, $\gamma$ could be negative even if $r$ is positive. Hence, the growth rate $\gamma$ is a better indicator of the long-term value of an investment than $r$.
Another way to see this is to notice that if $\gamma>0$, as $t$ increases the distribution of $S(t)$ becomes more skewed to the right with an ever greater right tail. As a consequence, as time passes the mean of the distribution moves deeper into the right tail, exponentially diverging from the center of the distribution, i.e. the median. To see this, notice that the ratio of the mean to the median is equal to $e^{\frac{1}{2}\sigma t}$.
Notice that the compounded rate of return $\gamma$ and the instantaneous rate of return $r$ $$\frac{d S(t)}{S(t)} = r t + \sigma dW_t$$ coincide for the value of a bank account with no risk paying interest rate $r$, i.e. $r = \gamma$ because in such case there is no randomness (aka risk) and hence $\sigma =0$. But for a risky investment, the compounded rate of return $\gamma < r$!
By applying Ito's Lemma on $f(t,x) = \ln x$ we get
$$d(\ln S_t) = \frac{1}{S_t}dS_t - \frac{1}{2}\frac{1}{S_t^2}dS_tdS_t$$
Substituting the expression for $dS_t$ we obtain
$$d(\ln S_t) = (\mu-\frac{1}{2}\sigma^2)dt + \sigma dW_t$$
This is an expression that you can integrate in a straightforward way if you know $\mu$ and $\sigma$. That leads to the expression for $S_t$ you had in your question.
Best Answer
To calculate the skewness (in fact, every order moments of GBM), we need Ito's lemma, https://en.wikipedia.org/wiki/It%C3%B4%27s_lemma, which states that for an Ito diffusion $$ dS_t =\mu_tdt+\sigma_tdW_t $$ and a twice-differentiable scalar function $F(t, S)$, we have $$ \begin{align} F(t,S_t)=F(0, S_0)&+\int ^t _0 F_1 (s, S)ds\\ &+\int ^t _0 F_2 (s, S)\mu_sds\\ &+\int ^t _0 F_2 (s, S)\sigma_sdW_s\\ &+\dfrac {1}{2}\int ^t _0 F_{22} (s, S)\sigma_s ^2ds\\ \end{align} $$ for all $t$. These four integrals come from Taylor expansion and the property of Brownian motion that $(dW_t)^2=dt$.
For geometric Brownian motion, we see that $$ \mu_t =\mu S_t, \sigma_t =\sigma S_t $$ moreover, to calculate skewness, we let $F(t, S)=S^3$, which doesn't depend on $t$. Then Ito's lemma gives $$ S^3_t=S_0^3 + \int_0^t (3\mu +3\sigma^2)S_s^3ds+ \int_0^t 3\sigma S_s^3dW_s $$ taking expectations, we get $$ \mathbb E S^3_t= S_0^3+ \int_0^t (3\mu +3\sigma^2)\mathbb E S_s^3ds $$ because the expectation of stochatic integral is always zero (provided $S_s^3$ is integrable).
Now let $h(t)=\mathbb E S^3_t$, then the above equation is $$ h(t)=S_0^3+(3\mu +3\sigma^2)\int_0^t h(s)ds $$ it's obvious: $$ h(t)=S_0^3 e^{(3\mu +3\sigma^2)}=\mathbb E S^3_t $$
Now we have all gradients for skewness: $\mathbb E S_t, \operatorname{Var} S_t, \mathbb E S^3_t$: $$ \mathbb E [S_t-\mathbb E S_t]^3=\mathbb E S_t^3-3\mathbb E S_t\operatorname{Var} S_t-(\mathbb E S_t)^3 $$ and the result follows: $$ S_0^3 e^{3\mu t}[e^{3\sigma^2 t}-3e^{\sigma^2 t}+2] $$