Let $\zeta(s) = \prod\limits_p (1 - p^{-s})^{-1}$ be the Riemann zeta function. If we define $L(s) = \pi^{-\frac{s}{2}}\Gamma(\frac{s}{2})\zeta(s)$, then $L$ admits a meromorphic continuation to the complex plane satisfying the functional equation
$$L(s) = L(1-s)$$
The proof techniques are well established, but at the time were very sophisticated, making use of the Fourier transform, Poisson summation, and complex line integrals.
My question is, how could anyone have known to thought to combine the gamma function and the Riemann zeta function in this way?
I see that the proof works, but I don't see how anyone could have ever thought of doing this. The same idea has carried over to prove the meromorphic continuation of other L-functions, but Riemann's combination of the gamma and zeta function seems to be the first of its kind.
Skimming the details of the proof since you seem to be familiar with it: the starting point to the classical proof of the functional equation is the expression $$\zeta(s) \Gamma(s) = \int_0^\infty \frac {t^{s-1}} {e^t-1} \, dt. \tag{$*$} $$ Already we can see the presence of the gamma function. Unfortunately this still only converges for $\text{Re}(s)>1$, so we haven't done any better than our original sum form. However, Riemann had the idea of deforming this into a contour integral $$\zeta(s) \Gamma(s) = \frac i {2 \sin (\pi s)} \int_C \frac {(-z)^{s-1}} {e^z-1} \, dz,$$ with the contour $C$ coming from positive infinity along the real line slightly in the upper half-plane, circling around the origin in the positive (counterclockwise) direction, and returning to positive infinity along the real line in the lower half-plane. This now looks like something we could evaluate using the residue theorem. In fact the contour as it is presents us with some challenges, but by deforming it suitably we can in fact do just this, and it turns out that the singularities will be at nonzero integral multiples of $2 \pi i$ with residues proportional (up to a function of $s$ independent of the integer $n$ in question) to $n^{1-s}$. Summing over these, we formally have a relationship between $\zeta(s) \Gamma(s)$ and $\zeta(1-s)$, which is excitingly close to what we're looking for; in fact filling in the factors we skipped over we get $$\zeta(1-s) = 2^{1-s} \pi^{-s} \cos \left( \frac {\pi s}2 \right) \zeta(s) \Gamma(s).$$ This is a perfectly fine form for a functional equation, but to get $L$ properly you can apply the duplication formula (I think Legendre's, definitely before Riemann). I was going to write out the derivation, but realized that a) it's not terribly relevant to the question and b) I'm bad at algebra, but it can be done.
As to where the initial equation $(*)$ comes from or why it's natural: if you start with the gamma function, by a change of variables you can naturally involve a factor of $n^{-s}$; and at this point to a mathematician of the nineteenth century the obvious thing to do is to sum over $n \ge 1$, giving us the formula.