In this answer, I will present you a proof which came to my mind five years ago when I was meeting people in a Waldorf session in my hometown...
First, let us write a permutation $\pi$ on $n$ numbers $\left[n\right]$ as
$$\pi = \begin{pmatrix}
1 & 2 & 3 & \cdots & n \\
\pi(1) & \pi(2) & \pi(3) & \cdots & \pi(n)
\end{pmatrix}$$
Next, we will define a function $C(\pi)$ as the number of cycles per given permutation $\pi$. The set of all permutations of the length $n$ is denoted as $S_n$. For length of a permutation we will write $l(\pi)=n$ Then, let us consider this sum
$$C_n(x) = \sum_{l(\pi)=n} x^{C(\pi)}$$
where we sum over all permutations $\pi$ with length $n$. Clearly $C_n(1)=n!$ and $C_1(x) = x$. We will show that there is a simple recurrence relation namely $C_{n+1}(x)=(x+n)C_n(x).$
First, let us construct the set of all permutations on $\left[n+1\right]$. First, choose any $\pi$ with $l(\pi)=n$ and add a one-cycle $(n+1, n+1)$, i.e. we have the permutation $\pi^*$
$$\pi^* = \begin{pmatrix}
1 & 2 & 3 & \cdots & n & n+1 \\
\pi(1) & \pi(2) & \pi(3) & \cdots & \pi(n) & n+1
\end{pmatrix}$$
This permutation has therefore exactly $C(\pi^*) = C(\pi)+1$ cycles. From this permutation, using one inversion, we can construct another $n$ distinct permutations ${\pi^*_i}_{i=1}^n$ of the form
$$\pi^*_i = \begin{pmatrix}
1 & 2 & 3 & \cdots & i & \cdots & n & n+1 \\
\pi(1) & \pi(2) & n+1 & \cdots & n & \cdots & \pi(n) & \pi(i)
\end{pmatrix}$$
The number of cycles of $\pi^*_i$ is one less than $\pi^*$ since the inversion merged two disjoint cycles, therefore $C(\pi^*_i)=C(\pi)$ (see the image below for an explanation, there $\pi^*$ stands for a composition of $\pi$ with the inversion of elements $x,y$).
We have not disrupted, however, the order of numbers ${1,2,3,\dots,n}$ so there exist an injection from $\pi^*_i$ to $\pi$. The total number of $\pi^*_i$ permutations is therefore $n n!$, together with $\pi^*$ permutations we have constructed $nn!+n! = (n+1)!$ distinct permutations so we have covered all permutations in $S_{n+1}$. Therefore
$$C_{n+1}(x) = \!\!\!\!\sum_{l(\pi)=n+1} \!\! \!\! x^{C(\pi)} = \!\!\!\!\sum_{l(\pi^*)=n+1} \!\! \!\!\! x^{C(\pi^*)} + \!\sum_{i=1}^{n}\!\sum_{l(\pi^*_i)=n+1} \!\! \!\!\! x^{C(\pi^*_i)} = \!\sum_{l(\pi)=n} \!\! x^{C(\pi)+1} + \!\sum_{i=1}^{n}\!\sum_{l(\pi)=n+1} \!\! \!\! x^{C(\pi)}= (x+n)C_n(x)$$
I.e., this shows that
$$C_n(x) = (x+n-1)(x+n-2)\cdots (x+3)(x+2)(x+1)x$$
or
$$D_n(x) = \ln C_n(x) = \sum_{k=1}^n \ln(x+k-1)$$
Let us define the nth Harmonic number of mth kind $H_n^{(m)}$ as
$$H_n^{(m)} = \sum_{k=1}^{n} \frac{1}{k^m}$$
with assymptotics $H_n=H_n^{(1)} \sim \gamma + \ln n$, where $\gamma$ is the Euler-Mascheroni constant, and for $m>1$ we have in the limit
$$\lim_{n\rightarrow\infty}H_n^{(m)} = \sum_{k=1}^{\infty} \frac{1}{k^m} = \zeta(m)$$
where $\zeta$ is the Riemann zeta function. Clearly, you can check for yourself, for derivatives of $D_n(x)$ at $x=1$, that
$$D_n^{(m)}(1) = (-1)^{m-1}(m-1)! H_n^{(m)}$$
for $m>1$ and $D_n(1)=\ln n!$. The complete statistic is given by $C_n(x)$ (or $D_n(x)$). To find the mean $\mu_n$, an average number of cycles in a random permutation of length $n$, we just have to compute the ratio
$$\mu_n=\mathbb{E}\left[C\right]=\frac{1}{n!}\sum_{l(\pi)=n} C(\pi) = \frac{1}{n!}\left(\sum_{l(\pi)=n} x^{C(\pi)}\right)'\bigg{|}_{x=1}=\frac{1}{n!}\left(C_n(x)\right)'=\frac{C_n(1)}{n!}D_n'(1)=H_n=\sum_{k=1}^n \frac{1}{k}.$$
The variance is defined $\sigma^2_n = \mathbb{E}\left[C^2\right]-\mathbb{E}^2\left[C\right]$. First we have
$$\mathbb{E}\left[C^2\right]=\frac{1}{n!}\sum_{l(\pi)=n} C^2(\pi) = \frac{1}{n!} (x C'_n(x))'\bigg{|}_{x=1}=\frac{1}{n!} (x C_n(x) D'_n(x))'\bigg{|}_{x=1}=\frac{1}{n!} C_n(1)\left(D'_n(1)+D'^2(1)+ D''_n(1)\right)= H_n^2 + H_n + H_n^{(2)}.$$
Hence,
$$\sigma_n^2 = H_n + H_n^{(2)} \rightarrow \gamma + \frac{\pi^2}{6} + \ln n$$
