How can we derive from $\{G,H\}=0$ that $G$ generates a transformations which leaves the form of Hamilton's equations unchanged? In the Hamiltonian formalism, a symmetry is defined as transformation generated by a function $G$ is a symmetry if
$$\{G,H\}=0 ,$$
where $H$ denotes the Hamiltonian. 
On the other hand, a symmetry is a transformations which map each solution of the equation of motion into another solution. And this requires that the form of the equation of motion remains unchanged. 
Therefore, it should be possible to show that it follows from $\{G,H\}=0 $ that Hamilton's equations are unchanged by the transformation generated by $G$. 
Concretely, we have
\begin{align}
q \to q' &= q +  \epsilon \frac{\partial G}{\partial p} \\
p \to p' &= p -   \frac{\partial G}{\partial q} \\
H \to H' &=H + \{H,G\} 
\end{align}
and we want to show that if for the original $q$ and $p$ Hamilton's equations 
\begin{align}
 \frac{dp}{dt}&=  -\frac{\partial H}{\partial q} \\
\frac{dq}{dt}  &=  \frac{\partial H}{\partial p} 
\end{align}
hold, they also hold for $q'$ and $p'$:
\begin{align}
 \frac{dp'}{dt}&=  -\frac{\partial H'}{\partial q'} \\
\frac{dq'}{dt}  &=  \frac{\partial H'}{\partial p'} 
\end{align}
How can this be shown explicitly?

Using the transformation rules explicitly yields for Hamilton's first equation
\begin{align}
 \frac{dp}{dt}&=  -\frac{\partial H}{\partial q} \\
\therefore \quad  \frac{d(p' +   \frac{\partial G}{\partial q})}{dt}&=  -\frac{\partial (H + \{H,G\} )}{\partial (q' +  \epsilon \frac{\partial G}{\partial q} )} \\
\therefore \quad  \frac{d(p' +   \frac{\partial G}{\partial q})}{dt}&=  -\frac{\partial H }{\partial (q' +  \epsilon \frac{\partial G}{\partial q} )} \\
\end{align}
But I've no idea how to proceed from here.
 A: A generating function$^1$ $\epsilon G(q,p,t)$ [where $\epsilon$ is an infinitesimal parameter] can to first order in $\epsilon$ be identified with a generating function $\epsilon G(q,P,t)$ for a type 2 canonical transformation (CT), $P=p+{\cal O}(\epsilon)$, cf. Ref. 1.
On the other hand, a CT takes Hamilton's equations into Kamilton's equations, with
$$K(Q,P,t)~=~H(q,p,t)+\epsilon \frac{\partial G(q,P,t)}{\partial t}, $$ 
so yes, it leaves Hamilton's equations form-invariant.
References:


*

*H. Goldstein, Classical Mechanics; eqs. (9.61)-(9.63).


--
$^1$ $G$ and $H$ do not have to Poisson commute. 
A: Express your equations of motion as 
$$
\dot q= [H,q]\\\\\dot p=[H,p]
$$
Note that, on the mass shell, any function $f(q,p)$ obeys $\dot f(q,p)=[H,f]$. Now just hit the commutator $[G, ]$ on both sides in the equations of motion. Since $[G, H]=0$ then $G$ commute with the time derivative, i.e. $\dot G=0$. Using the Jacobi identity on the right hand side gives you 
$$
[G,[H,q]] = [H,[G,q]] + [[G,H],q] = [H,[G,q]]
$$
Now you are finish since you can use the linearity of $[H,]$ to add this new equation into your old equation
$$
\frac{d}{dt}(q+\epsilon [G,q])=[H,(q+\epsilon [G,q])]
$$
The same is true for the $p$-equation.
A: You just might have gotten your logic crossed. A generic canonical transformation, in your language, is generated by a $G(q,p)$ which need not PB-commute with the hamiltonian, and so need not be a symmetry (even though it might, if you felt like it.).
$$
q'=q+\epsilon \frac{\partial G}{\partial p}+...= q-\epsilon \{ G,q \}+ ... = e^{-\epsilon \{ G , ~\cdot        }   ~q ,\\
p'=p-\epsilon \frac{\partial G}{\partial q}+...= p-\epsilon \{ G,p \}+ ... = e^{-\epsilon \{ G , ~\cdot        }   ~p ,  \\
H(q'p')= H(q,p) -\epsilon \{ G, H \} +... ,\\
G(q',p')= G(q,p).
$$
You may thence confirm that $\{ q'  , p' \}=1$, and vice versa, i.e. transforming backwhards from primed to unprimed.  In addition, $H(q', p')=H(q,p)$ if the Generator PB-commutes with the Hamiltonian, a special case indeed. 
Start by thinking of the primed variables as functions of the unprimed ones. To get the unprimed evolution equations of primed quantities qua arbitrary functions , 
$$
\dot q'= - \{ H,q' \}  = - \{ H(q',p') +\epsilon \{ G,H\} ,q' \} + O(\epsilon^2) \\
\dot p '=-\{ H, p'\}  = - \{ H(q',p') +\epsilon \{ G,H\} ,p' \} + O(\epsilon^2) ,
$$
where, to this order in $\epsilon$, it does not matter which variables (primed or unprimed) we use in the H in the term multiplied by $\epsilon$ (!), and of course it never mattered which variables we use in G. Moreover, it does not matter what variables we use for the PBs, as these preserve the PBs: the magic of the unit Jacobian (Liouville's theorem). 
Now elect to consider only the primed variables, observing that, for $K(q',p')\equiv H(q,p)$ $= H(q',p')+\epsilon \{G, H  \} $, 
Hamilton's equations hold for everything, including PBs, in  primed variables,
$$
\dot q'= - \{ K,q' \} ' \qquad 
\dot p '=-\{ K, p'\} ' .
$$
In the very special case when $\{ G,H   \} =0$, you see that $K(x,y)=H(x,y)$, the hamiltonian is invariant. 
I can't urge you enough to  examine a toy example such as $2G=p^2+q^2$, work out the all-orders canonical transformation 
$$
q'=\cos \epsilon  ~ q   +\sin \epsilon ~p ,\\
p'= \cos \epsilon ~p   -\sin \epsilon ~q,
$$
and consider the dynamics of a hamiltonian for which this is a symmetry, $ G$, and one where it is not, $(p^2+q^4)/2$.
