The Riemann zeta function is given for $Re s>1$
$$\zeta(s) = \prod\limits_{p \; prime} ( 1-p^{-s}) $$
This product converges absolutely in $Re s >1$, hence it does not vanish in $Re s>1$. Actually the product also converges locally uniformly, which implies that $\zeta(s)$ is holomorphic for $Re s>1$.
The functional equation follows from the Poisson summation formula, which is a Fourier theoretic argument. The functional equation is given by $$ \Lambda(s) = \pi^{-s/2} \Gamma(s/2) \zeta(s) = \Lambda(1-s).$$ This implies that all zeros lie inside $0 \leq Re s \leq 1$. We will need to evaluate a certain integral $$ * = \int\limits_{Re\; s = 1 + \epsilon/2 \cup Re \; s = -\epsilon/2} G(s) \frac{\Lambda'}{\Lambda} (s) d s$$
We have now for a meromorphic function $F$ and a holomorphic function $G$, that the contour integral $$\frac{1}{2 \pi i} \int\limits_{C} G(z) \frac{F'}{F}(z) d z = \sum\limits_{\rho \; zero \; of \; F} G( \rho) - \sum\limits_{\nu \; pol} G( \nu).$$ This can be derived in the same lines as the argument principle,and does not need Hadamard factorization theorem.
Apply this to the function $F = \Lambda$ and $G$ being holomorphic in $ - \epsilon \leq Re\; s \leq 1$ with $|G(z)| \ll (1+|z|^2)^{-1-\epsilon'}$ for some $\epsilon, \epsilon' >0$. We choose the contour $C=C(T)$ being the boundary of $ - \epsilon/2 < Re\;s < 1 + \epsilon/2$ and $| Im \;s | \leq T$, where \zeta does not vanish on $Im \; s = \pm T$. This give an expression for $*$ involving the nontrivial zeros of $\zeta$.
Using the Euler product, we can also derive a nice explicit expression $$ \frac{\Lambda'}{\Lambda} (s) = -1/2 \log \pi -\frac{1}{2} \frac{\Gamma'}{\Gamma}(s/2) + \sum\limits_{p \; prime} \frac{\log p}{p^{-s}}.$$ This gives an expression for $*$ involving the primes.
Assuming certain boundedness conditions on $\Lambda(s)$ in $0 \leq Re \; s \leq 1$, we are allowed to choose $C(T)$ with $T \rightarrow \infty$ and derive the formula. The boundedness conditions follow from the Hadamard three lines principle or the Phragmen Lindeloeff principle (this is not entire function theory, but this only a complex analysis argument), then the explicit formula follows.
You are right that the entire function theory implies that $\Lambda$ has necessary many zeros, since it is of exponential type $1$ because of the factor $\Gamma$. If you want to derive this without using the merophorphicity of $\Lambda$, you might want to try to deduce this without knowledge over the primes and by inserting an approppiate chosen test function in the explicit formula. I have never seen this been worked out, but choosing an appropiate function $g$ being supported in $ - \log 2 < t \log 2$ etc. should lead to a rough asymptotic of the zeros without any information used about the primes, but possibly a weaker error term. Look at similiar techniques used in Werner Mueller and Erez Lapid's article Chapter 2 of http://www.math.uni-bonn.de/people/mueller/papers/orbint09.11.pdf for the Weyl law. The Selberg trace formula has many analogies with the explicit formula.
One interesting, but technical derivation of the explicit formula using only the languages of the adeles, harmonic analysis and no entire function theory at all was given by Ralf Meyer: http://arxiv.org/pdf/math/0311468v3
However, the Fourier transform of function with certain growth properties have some holomorphicity conditions, so you basically will just hide the complex analysis arguments, but be able to deduce the results from Fourier analysis only.