The first step is to recognize the fraction involving $e^{\pi\sqrt{n^2+m^2}}$ as $\coth (\pi\sqrt{n^2+m^2})$. Then we can use the Mittag-Leffler expansion of $\coth$, $$ \pi\coth(\pi z) = \frac{1}{z} + 2 \sum_{n=1}^\infty \frac{z}{z^2+n^2}\ . $$ Applying this to $z = \sqrt{m^2+n^2}$, the stuff inside the parantheses in $I$ can be rewritten as $$ \frac{1}{\sqrt{m^2+n^2}} + 2\sum_{p=1}^\infty \frac{\sqrt{m^2+n^2}}{m^2+n^2+p^2} - \frac{1}{\sqrt{m^2+n^2}}\ . $$ After canceling out the $1/\sqrt{m^2+n^2}$, note that the $\sqrt{m^2+n^2}$ factor outside of the parentheses also gets canceled. Thus $$ \frac{\pi^6}{4}I^{-1} = \sum_{m=1}^\infty \sum_{n=1}^\infty \sum_{p=1}^\infty \frac{2}{m^2n^2(m^2+n^2+p^2)}\ . $$ To evaluate this sum, note that if we permute the dummy variables $m\mapsto n$, $n\mapsto p$, $p\mapsto m$, the summand can be rewritten as $$ \frac{2}{n^2p^2(m^2+n^2+p^2)}\ . $$ Permutting again turns the summand into $$ \frac{2}{p^2m^2(m^2+n^2+p^2)}\ . $$ Since relabeling the dummy variables clearly doesn't change the value of the series, we could add all three representations together and get $$ \begin{align} 3\left(\frac{\pi^6}{8I}\right) &= \sum_{m=1}^\infty \sum_{n=1}^\infty \sum_{p=1}^\infty \left(\frac{1}{m^2n^2(m^2+n^2+p^2)} +\right.\\ & \quad\left. \frac{1}{n^2p^2(m^2+n^2+p^2)} + \frac{1}{p^2m^2(m^2+n^2+p^2)}\right) \\ &= \sum_{m=1}^\infty \sum_{n=1}^\infty \sum_{p=1}^\infty \frac{p^2+m^2+n^2}{m^2n^2p^2(m^2+n^2+p^2)}\\ & = \sum_{m=1}^\infty \sum_{n=1}^\infty \sum_{p=1}^\infty \frac{1}{m^2n^2p^2} \\ & = \left(\sum_{m=1}^\infty \frac{1}{m^2}\right)^3 \\ & = \left(\frac{\pi^2}{6}\right)^3 \end{align} $$ From this we can easily solve for $I = 81$.