As indicated by Igor Rivin, the volume of the unitary group is given by 
$vol(U(N))=(2\pi)^{(N^2+N)/2}/\prod_{k=1}^{N-1} k!$.

The denominator is the Barnes G-function, which is well-known : 

http://en.wikipedia.org/wiki/Barnes_G-function

and in particular has a known Stirling-like asymptotic expansion for large $N$:

$log(\prod_{k=1}^{N-1} k!) \sim$
$\frac{N^2}{2} log N - \frac{1}{12} log N - \frac{3}{4}N^2
+\frac{N}{2} log(2 \pi)
+ \zeta'(-1) + \sum_{g \geq 2} \frac{B_{2g}}{2g(2g-2)} N^{2-2g}$.

Comparing with the Harer-Zagier formula $\chi(M_g)=\frac{B_{2g}}{2g(2g-2)}$, we obtain
$vol(U(N))\sim$
$ - \frac{N^2}{2} log(N) + \frac{N^2}{2}(log(2 \pi)+\frac{3}{2}) + \frac{1}{12}
log(N) - \zeta'(-1)-\sum_{g \geq 2} \chi(M_g) N^{2-2g}$

which, up to probable typos and forgeted terms, is the asymptotic expansion of the question.

Of course, in such a proof, the fact that unitary groups and moduli spaces of Riemann surfaces are related appears as a coincidence: essentially we have just taken two places in mathematics where Bernoulli numbers appear. We can ask if there is a more direct intrinsic explanation of this relation. I don't think that such explanation is known at the level of rigorous mathematics but one is known at the level of theoretical physics. On general grounds, it is expected that gauge theories of group $U(N)$ are related in the large $N$ limit to a form of string theory. The first argument in this direction was given by 't Hooft in the 70's and is the observation that Feynman diagrams in perturbative $U(N)$ gauge theory can be rewritten as double-line graphs, or ribbon graphs, that it is possible to obtain closed surfaces from ribbon graphs by gluing disks along their boundary components, and that in some appropriate limit the series of Feynman diagrams organizes as a genus expansion of these surfaces.

In fact, it is possible to prove the Harer-Zagier relation along these lines by describing the moduli space of Riemann surfaces in terms of ribbon graphs, interpreting these ribbon graphs as the perturbative expansion of some $N$ by $N$ matrix model, solving this matrix model, which gives something containing the $\Gamma$ function and then expanding the solution. In this proof, which can be found in an appendix to Kontsevich's paper on Witten's conjecture, http://www.ihes.fr/~maxim/TEXTS/intersection_theory_6.pdf , the Bernoulli numbers appearing in the Harer-Zagier formula really comes from the Stirling expansion of the $\Gamma$-function.

Making 't Hooft idea concrete is one of the main theme of modern string theory and can go under various more or less general and more or less precise names: gauge/gravity duality, AdS/CFT correspondence, open/closed duality, holographic relation... One explicit example of that is the Gopakumar-Vafa correspondence asserting the equivalence of Chern-Simons theory of group $U(N)$ and level $k$ on the 3-sphere with the A-model of the topological string, i.e. Gromow-Witten theory, on the resolved conifold, i.e. the total space of $\mathcal{O(-1)}\oplus {\mathcal{O(-1)}}$ over $\mathbb{P}^1$, with $\mathbb{P}^1$ being of volume $t=\frac{2 \pi N}{k+N}$ and with a string coupling constant 
$g_s = \frac{2 \pi}{k+N}$. As the volume of $U(N)$ appears explicitely in the one-loop perturbative expansion of Chern-Simons theory on the 3-sphere, it is possible to "explain" the asymptotics expansion of these volumes in terms of moduli spaces of Riemann surfaces. 
Of course, all that is not a proof and the direct matching of the two sides of the equalities is often used as a support of physicists conjectures but I wanted to mention it because it is a natural circle of ideas in which the formula of the question appears naturally.