We first notice that
\begin{align*}
\int_{0}^{\infty}\frac{\operatorname{gd}(x)}{e^x-1}\,dx
&= \int_{0}^{\infty} \frac{1}{e^x-1} \left( \int_{0}^{x} \frac{dy}{\cosh y} \right) \, dx \\
&= \int_{0}^{\infty} \frac{1}{\cosh y} \left( \int_{y}^{\infty} \frac{dx}{e^x - 1} \right) \,d y \\
&= -2 \int_{0}^{\infty} \frac{\log(1 - e^{-y})}{e^y + e^{-y}} \, dy \\
&=-2\int_{0}^{\frac{\pi}{4}}\log(1-\tan\theta)\,d\theta \tag{$e^{-x}=\tan\theta$}.
\end{align*}
The last integral is our starting point. We introduce two tricks to evaluate this.
Step 1. Notice that $\tan(\frac{\pi}{4}-\theta)=\frac{1-\tan\theta}{1+\tan\theta}$. So by the substitution $\theta \mapsto \frac{\pi}{4}-\theta$, it follows that
$$ \int_{0}^{\frac{\pi}{4}}\log(1+\tan\theta)\,d\theta
= \int_{0}^{\frac{\pi}{4}}\log\left(\frac{2}{1+\tan\theta}\right)\,d\theta $$
and hence both integrals have the common value $\frac{\pi}{8}\log 2$. Applying the same idea to our integral, it then follows that
\begin{align*}
-2\int_{0}^{\frac{\pi}{4}}\log(1-\tan\theta)\,d\theta
&= -2\int_{0}^{\frac{\pi}{4}}\log\left(\frac{2\tan\theta}{1+\tan\theta}\right)\,d\theta \\
&= -2\int_{0}^{\frac{\pi}{4}}\log\tan\theta \, d\theta - \frac{\pi}{4}\log 2.
\end{align*}
Step 2. In order to compute the last integral, we notice that for $\theta\in\mathbb{R}$ with $\cos\theta\neq0$, we have
\begin{align*}
-\log\left|\tan\theta\right|
&= \log\left|\frac{1+e^{2i\theta}}{1-e^{2i\theta}}\right|
= \operatorname{Re} \log\left(\frac{1+e^{2i\theta}}{1-e^{2i\theta}}\right) \\
&= \operatorname{Re}\left( \sum_{n=1}^{\infty} \frac{1+(-1)^n}{n} e^{2in\theta} \right) \\
&= \sum_{k=0}^{\infty} \frac{2}{2k+1}\cos(4k+2)\theta.
\end{align*}
So by term-wise integration, we obtain
\begin{align*}
-2\int_{0}^{\frac{\pi}{4}}\log\tan\theta \, d\theta
&= \sum_{k=0}^{\infty} \frac{4}{2k+1} \int_{0}^{\frac{\pi}{4}} \cos(4k+2)\theta \, d\theta \\
&= 2 \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)^2}
= 2K,
\end{align*}
where $K$ is the Catalan's constant.
Best Answer
You may refer to the entry 9.10.17 of DLMF, which describes the Mellin transform of the Airy function $\operatorname{Ai}$:
$$ \int_{0}^{\infty} t^{\alpha-1}\operatorname{Ai}(t)\,\mathrm{d}t = \frac{\Gamma(\alpha)}{3^{(\alpha+2)/3}\Gamma\big(\frac{\alpha+2}{3}\big)}. \tag{9.10.17}$$
Differentiating both sides with respect to $\alpha$ 3-times and plugging $\alpha = 1$ gives the answer to your integral in terms of polygamma functions with further simplifications available.
Sketch of proof of $\text{(9.10.17)}$. We begin with the integral representation
$$ \operatorname{Ai}(x) = \frac{1}{\pi} \int_{0}^{\infty} \cos\left(\frac{t^3}{3} + xt \right) \, \mathrm{d}t. $$
Taking Mellin transform and switching the order of integration,
\begin{align*} \int_{0}^{\infty} x^{\alpha-1}\operatorname{Ai}(x)\,\mathrm{d}x &= \frac{1}{\pi} \operatorname{Re}\bigg[ \int_{0}^{\infty}\int_{0}^{\infty} x^{\alpha-1} e^{\frac{1}{3}it^3 + ixt} \, \mathrm{d}x\mathrm{d}t \bigg] \\ &= \frac{1}{\pi} \operatorname{Re}\bigg[ \int_{0}^{\infty} \frac{\Gamma(\alpha)}{(-it)^{\alpha}} e^{\frac{1}{3}it^3} \, \mathrm{d}t \bigg] \\ &= \frac{1}{\pi} \operatorname{Re}\bigg[ \Gamma(\alpha) i^{\alpha} \frac{\Gamma\big(\frac{1-\alpha}{3}\big)}{3\cdot(-i/3)^{\frac{1-\alpha}{3}}} \bigg] \\ &= \frac{1}{\pi} \operatorname{Re}\bigg[ \Gamma(\alpha) i^{\frac{2\alpha+1}{3}} \frac{\Gamma\big(\frac{1-\alpha}{3}\big)}{3^{\frac{\alpha+2}{3}}} \bigg] \\ &= \frac{\Gamma(\alpha)\Gamma\big(\frac{1-\alpha}{3}\big)}{\pi 3^{\frac{\alpha+2}{3}}}\cos\big( \tfrac{1}{6}\pi + \tfrac{1}{3}\alpha\pi \big). \end{align*}
Now applying the Euler's reflection formula $\Gamma(s)\Gamma(1-s) = \pi \csc(\pi s)$ with $s = \frac{1-\alpha}{3}$ yields the right-hand side of $\text{(9.10.17)}$.