How can this integral expression for the difference between two $\zeta(s)$s be explained?

With $z=\sigma + x \, i;\,\, \sigma,x \in \mathbb{R}$, numerical evidence strongly suggests that:

$$\displaystyle \int_{0}^{\infty} \,\zeta(z)\,\zeta(\overline{z})\,\Gamma(z)\,\Gamma(\overline{z}) \,dx =\pi\,\Gamma(2\,\sigma)\,\big(\zeta(2\,\sigma-1)-\zeta(2\,\sigma)\big)$$

for all $\sigma > 1$.

Could this be proven?

Just to share that the equation can be extended towards $0<\sigma<1$ by starting from: $$\Gamma(s)\, \zeta(s) = \int_0^\infty x^{s-1}\left(\frac{1}{e^{x}-1}-\frac{1}{x}\right)\,dx$$ and following the same logic as in the answer below. This gives the closed form:

$$\displaystyle \int_{0}^{\infty} \,\zeta(z)\,\zeta(\overline{z})\,\Gamma(z)\,\Gamma(\overline{z}) \,dx =\pi\,\Gamma(2\,\sigma)\,\left(\frac{(2\,\sigma-3)}{(2\,\sigma-1)} \,\zeta(2\,\sigma-1)-\zeta(2\,\sigma)\right)$$

Note that since $\sigma$ could also be $\in \mathbb{C}$, this implies that when $\rho=2\,\sigma$ or $\rho=2\,\sigma-1$ ($\rho$ = a non-trivial zero of $\zeta(s)$), then the RHS becomes fully multiplicative.

The process can be extended indefinitely by for instance starting from $$\Gamma(s)\, \zeta(s) = \int_0^\infty x^{s-1}\left(\frac{1}{e^{x}-1}-\frac{1}{x} +\frac12\right)\,dx$$ which gives for the domain $-1<\sigma<0$ the following closed form:

$$\displaystyle \int_{0}^{\infty} \,\zeta(z)\,\zeta(\overline{z})\,\Gamma(z)\,\Gamma(\overline{z}) \,dx =\pi\,\left(2\,\sigma-3)\,\Gamma(2\,\sigma-1)\,\zeta(2\,\sigma-1\right)$$

And subsequently for the domain $-2<\sigma<-1$ we start from: $$\Gamma(s)\, \zeta(s) = \int_0^\infty x^{s-1}\left(\frac{1}{e^{x}-1}-\frac{1}{x} +\frac12-\frac{x}{12}\right)\,dx$$

which gives the closed form:

$$\displaystyle \int_{0}^{\infty} \,\zeta(z)\,\zeta(\overline{z})\,\Gamma(z)\,\Gamma(\overline{z}) \,dx =\pi\,\Gamma(2\,\sigma)\,\left(\frac{(2\,\sigma-3)}{(2\,\sigma-1)} \,\zeta(2\,\sigma-1)-\frac{\sigma}{3}\,\zeta(2\,\sigma+1)\right)$$

Don't think there is an easy pattern or generic formula since Bernoulli numbers are involved.

Again, you have to look at the Fourier/Laplace transform. For $Re(s) > 1$ : $$\Gamma(s) \zeta(s) = \int_0^\infty \frac{x^{s-1}}{e^{x}-1}dx \overset{(x \, = \,e^u)}= \int_{-\infty}^\infty \frac{e^{su}}{e^{e^u}-1}du = \int_{-\infty}^\infty \frac{e^{-su}}{e^{e^{-u}}-1}du$$

i.e. $\Gamma(s)\zeta(s)$ is the Laplace transform of $\frac{1}{e^{e^{-u}}-1}$ and $F_\sigma(\xi) = \Gamma(\sigma+2i \pi \xi)\zeta(\sigma+2i \pi \xi)$ is the Fourier transform of $f_\sigma(u) = \frac{e^{-\sigma u}}{e^{e^{-u}}-1}$.

Now we can apply the Parseval theorem and get $$\int_{-\infty}^\infty |F_\sigma(\xi)|^2d\xi = \int_{-\infty}^\infty |f_\sigma(u)|^2du = \int_{-\infty}^\infty \frac{e^{-2\sigma u}}{(e^{e^{-u}}-1)^2}du $$

$$=\int_{-\infty}^\infty \frac{e^{2\sigma u}}{(e^{e^{u}}-1)^2}du = \int_0^\infty \frac{x^{2\sigma-1}}{(e^{x}-1)^2}dx = \int_0^\infty x^{2\sigma-1}\sum_{n=2}^\infty (n-1) e^{-nx}dx = (\zeta(2\sigma-1)-\zeta(2\sigma))\Gamma(2\sigma)$$

And note that $\zeta(s)$ is the Laplace transform of the distribution $\sum_{n=1}^\infty \delta(u-\ln n)$, while (under the Riemann hypothesis) the non-trivial zeros define the spectrum of the distribution $\sum_{p^k} \frac{\ln p}{\sqrt{p^k}} \ (\delta(u-\ln p^k)+\delta(u+\ln p^k))$