Rather than the problem of why the zeta function has non-trivial zeros, let me address
Gowers's question of why the error term in the prime number theorem needs to be large. The short answer that I propose is: because the integers are so well distributed. To make this precise, I shall prove a general result on semigroups, showing that
either the "integers" in the semigroup or the "primes" must be poorly distributed -- one may think of this as an "uncertainty principle" that both primes and integers cannot be simultaneously smoothly behaved. This has the flavor of Beurling generalized primes, but I don't recall this result in the literature; maybe it exists already (indeed it does, see the edit below). Also note that the proof will not make any use
of the functional equation for $\zeta(s)$ as this does not exist in a general semigroup.
EDIT: Indeed doing a literature search a few hours after posting this,
I found a paper of Hilberdink http://www.sciencedirect.com/science/article/pii/S0022314X04002069 which proves the
Theorem below, and with a similar method of proof.
Suppose that $1<p_1< p_2 <\ldots$ is a sequence of real numbers (the "primes"), and
$1=a_1 \le a_2 \le \ldots$ is the semigroup generated by them (the "integers"). Initially I made the assumption that the "integers" are distinct, and satisfy a mild spacing
condition $a_{n+1}-a_n \gg n^{-1}$ (so that in particular there is unique factorization), but this is not necessary. Let $N(x)$ denote the number of
"integers" below $x$, and
$$
P(x) = \sum_{p_n^{k} \le x} \log p_n,
$$
which is the analog of the usual $\psi(x)$ (counting prime powers with weight $\log p$). Assume that for some $\delta>0$
$$
N(x) = Ax +O(x^{\frac 12-\delta}),
$$
for some non-zero constant $A$, and that
$$
P(x) = x + O(x^{\frac 12-\delta}).
$$
Theorem. Either the asymptotic formula for $N(x)$ or the asymptotic
formula for $P(x)$ must fail.
Put
$$
\zeta_A(s) = \sum_{n=1}^{\infty} a_n^{-s} = \prod_{n} \Big(1-\frac{1}{p_n^s}\Big)^{-1}.
$$
By our assumptions on $N$ and $P$, the sum and product above converge absolutely in Re$(s)>1$. By the assumption on $N(x)$, $\zeta_A(s)$ extends to an analytic function in Re$(s)>1/2-\delta$ except for a simple pole at $s=1$ with residue $A$. By the assumption on $P$, we see that the logarithmic derivative
$$
-\frac{\zeta_A^{\prime}}{\zeta_A}(s) = \sum_{n, k} \frac{\log p_n}{p_n^{ks}}
$$
extends analytically to Re$(s)>1/2-\delta$ except for a simple pole at $1$. Thus
$\zeta_A(s)$ has no zeros in Re$(s)>1/2-\delta$.
New edit: Sketch of a second proof. Adapting the argument that the Riemann hypothesis implies the Lindelof hypothesis (see below), we obtain that $|\zeta_A(s)| \ll (1+|s|)^{\epsilon}$ provided $s$ is not close to the pole at $1$, and that Re$(s)>1/2-\delta/2$. From this and a standard contour shift argument we find that for large $N$ and any $t$,
$$
\sum_{a_n\le N} a_n^{it} = \frac{N^{1+it}}{1+it} +O(N^{1/2-\delta+\epsilon} (1+|t|)^{\epsilon}).
$$
What is used here is that we have Lindelof even a little to the left of the half line.
But the above identity can be seen to contradict the Plancherel formula. More precisely, let $T$ be a large power of $N$, and let $\Phi$ be a non-negative function supported in $[-1,1]$ with non-negative Fourier transform. Then we see that (discarding all but the diagonal terms)
$$
\int_{-\infty}^{\infty} \Big| \sum_{a_n\le N} a_n^{it}\Big|^2 \Phi(t/T) dt
\ge T{\hat \Phi}(0) \sum_{a_n\le N} 1 \sim TN {\hat \Phi}(0).
$$
On the other hand, if we use our identity then the above is seen to be
$$
\ll N^2 + T^{1+\epsilon} N^{1-2\delta}.
$$
This is a contradiction.
Original Proof: Below let's assume always that we are in the region Re$(s)>1/2-\delta$, and that
the imaginary part is large so that we are not near the pole at $1$. From
the analytic continuation of $\zeta_A$ (using that $N(x)$ is very regular), it follows
that there is an a priori polynomial bound $|\zeta_A(s)|\ll |s|^{B}$ in the region Re$(s)>1/2-\delta/2$. Thus there is a bound for the real part of $\log \zeta_{A}(s)$, and by the Borel-Caratheodory lemma (standard complex analysis) one can bootstrap this to a bound for $|\log \zeta_A(s)|$. Then applying the Hadamard three circle theorem to $\log \zeta_A(s)$ one obtains a much better bound: $|\zeta_A(s)| \ll |s|^{\epsilon}$. This is the usual proof that Riemann implies Lindelof. (At this stage, if we knew a "functional equation" we'd be done, as the usual $\zeta(s)$ is large when Re$(s)<1/2$. This point appeared in Matt Young's answers earlier.)
Knowing the Lindelof hypothesis for $\zeta_A(s)$ in Re$(s)> 1/2-\delta/2$, we can
show the following approximate formula: for $\sigma > 1/2- \delta/4$
$$
\zeta_A(\sigma +it) = \sum_{a_n \le N} a_n^{-\sigma-it} + O(|t|^{\epsilon} N^{-\delta/4}).
$$
The proof is standard; see the penultimate chapter of Titchmarsh for the real $\zeta(s)$ where this holds when $\sigma$ is strictly bigger than $1/2$, and our stronger result is true because we have Lindelof in a wider region.
Now we are ready to get our contradiction. Consider for large $T$
$$
\int_T^{2T} |\zeta_A(1/2+it)|^2 dt.
$$
To do this carefully it may be helpful to put in a smooth weight $\Phi(t/T)$ above
(but this is not a paper!). Using the approximate formula derived above, and our mild spacing condition $a_{n+1}-a_n \gg n^{-1}$,
we may see that for any $T^{\epsilon} \le N \le T^{1/10}$ we have
$$
\int_T^{2T} |\zeta_A(1/2+it)|^2 dt \sim T \sum_{a_n \le N} a_n^{-1} \sim AT \log N.
$$
But that's absurd! This completes our sketch proof.