Comenzando con
$$S(x) = \sum_{n\ge 1}\sum_{m\ge 1} e^{-mnx}$$
tenemos que para $x>0$ $m,n\ge 1$ el plazo $e^{-mnx} < 1$, por lo que podemos simplificar el interior de la suma para obtener
$$S(x) = \sum_{n\ge 1} \frac{e^{-nx}}{1-e^{-nx}}.$$
La suma plazo es armónico y puede ser evaluado por la inversión de sus Mellin transformar.
Recordar que el armónico suma de identidad
$$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$
donde $g^*(s)$ es la transformada de Mellin $g(x).$
En el presente caso tenemos
$$\lambda_k = 1, \quad \mu_k = k \quad \text{y} \quad
g(x) = \frac{e^{-x}}{1-e^{-x}}.$$
Necesitamos la Mellin transformar $g^*(s)$ $g(x)$ que es
$$g^*(s) = \int_0^\infty \frac{e^{-x}}{1-e^{-x}} x^{m-1} dx
= \int_0^\infty \sum_{q\ge 1} e^{-qx} x^{m-1} dx
\\= \sum_{q\ge 1} \int_0^\infty e^{-qx} x^{m-1} dx
= \Gamma(s) \sum_{q\ge 1} \frac{1}{q^s} = \Gamma(s) \zeta(s).$$
De ello se desprende que la Mellin transformar $Q(s)$ de la suma de armónicos $S(x)$ está dado por
$$Q(s) = \Gamma(s)\zeta(s)^2
\quad\text{porque}\quad
\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \zeta(s).$$
El Mellin de inversión integral aquí es
$$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$
de la cual evaluamos desplazando a la izquierda por una expansión de alrededor de cero.
La contribución de la pole en $s=1$ es
$$\mathrm{Res}\left(Q(s)/x^s; s=1\right) = - \frac{\log x}{x} + \frac{\gamma}{x}.$$
El polo en $s=0$ contribuye
$$\mathrm{Res}\left(Q(s)/x^s; s=1\right) = \frac{1}{4}.$$
Un simple cálculo muestra que el resto de los postes de contribuir
$$\sum_{q\ge 1} \frac{(-1)^q}{q!} \frac{B_{q+1}^2}{(q+1)^2} x^q.$$
Esta fórmula fue escogido por su simplicidad, otra fórmula puede ser obtenido mediante la observación de que los números de Bernoulli se desvanecen en números enteros impares.
La conclusión es que en una vecindad de cero, tenemos
$$S(x) \sim - \frac{\log x}{x} + \frac{\gamma}{x} + \frac{1}{4}
+ \sum_{q\ge 1} \frac{(-1)^p}{q!} \frac{B_{p+1}^2}{(p+1)^2} x^p$$
que es
$$S(x) \sim - \frac{\log x}{x} + \frac{\gamma}{x} + \frac{1}{4}
-{\frac {1}{144}}\,x-{\frac {1}{86400}}\,{x}^{3}-{\frac {1}{7620480}}\,{x}^{5}-
{\frac {1}{290304000}}\,{x}^{7}
-\cdots.$$
Como consecuencia de esta expansión, tenemos que
$$\lim_{x\to 0} \left(S(x)+\frac{\log x}{x}\right)\times x = \gamma,$$
lo que demuestra la conjetura.
Adenda. El residuo de $Q(s)$ en el poste de $s=1$ se calcula a partir de las tres de la serie de expansiones
acerca de $s=1$ $$\zeta(s)^2 = \frac{1}{(s-1)^2} + \frac{2\gamma}{(s-1)}+\cdots,
\quad \Gamma(s) = 1 - \gamma (s-1)+\cdots$$
y $$1/x^s = \frac{1}{x}-\frac{\log x}{x} (s-1)+\cdots.$$
