Before getting started, it's worth mentioning that you don't actually need twelve constants of motion, you only need six, one per degree of freedom. And you get that total energy is conserved for free, so you actually only need five.
Anyway, we know this is the Hamiltonian for the two body problem:
$$H = \frac{m_1}{2}(\dot x_1^2 +\dot y_1^2 +\dot z_1^2) +\frac{m_2}{2}(\dot x_2^2+ \dot y_2^2+ \dot z_2^2) + \frac{Gm_1m_2}{\sqrt{(x_1-x_2)^2+(y_1 - y_2)^2+(z_1-z_2)^2}}$$.
To start with, we'll cut six degrees of freedom down to three by switching to Jacobi coordinates:
$$\mathbf{r} = \mathbf{x_1}-\mathbf{x_2} \quad and\quad \mathbf{R} = \frac{m_1}{m_1+m_2}\mathbf{x_1}+\frac{m_2}{m_1+m_2}\mathbf{x_2}$$
Where $\mathbf{x_1}$ and $\mathbf{x_2}$ are the bodies' positions in Cartesian coordinates, $\mathbf{R}$ is the location of the center of mass, and $\mathbf{r}$ is the separation between the two bodies. With some quick algebra, the Hamiltonian turns into:
$$H = \frac{1}{2}(\dot{\mathbf{R}}^2 + \frac{m_1m_2}{m_1+m_2}\dot{\mathbf{r}}^2) + \frac{Gm_1m_2}{|\mathbf{r}|} $$
Hamilton's equations for $\mathbf{R}$ are
$$\frac{\partial H}{\partial \mathbf{p}_{\mathbf{R}}} = \dot{\mathbf{R}}\quad\rightarrow \quad\mathbf{p}_{\mathbf{R}} = \dot{\mathbf{R}} $$
and
$$\frac{\partial H}{\partial \mathbf{R}} = \dot{\mathbf{p}}_{\mathbf{R}} = 0\quad\rightarrow \quad\ddot{\mathbf{R}} = 0\quad$$
i.e. the velocity of the center of mass is constant. This accounts for three conserved quantities. The dynamics are the same if H varies by a constant, so we can drop the $\frac{1}{2}\dot{\mathbf{R}}^2$ term and we just have
$$H = \frac{m_1m_2}{2(m_1+m_2)}\dot{\mathbf{r}}^2 + \frac{Gm_1m_2}{||\mathbf{r}||} $$
Which is basically the Hamiltonian for a single particle in a stationary central potential, so we've succeeded in eliminating three degrees of freedom. Notice that $||\mathbf{r}||$ and $\dot{\mathbf{r}}^2 = ||\dot{\mathbf{r}}||^2$ are both invariant under a rotation in three dimensions, so the Hamiltonian (and therefore the Lagrangian) is invariant under such a transformation as well. A rotation in three dimensions can be decomposed into three rotations - one around the $x$ axis, one around the $y$ axis, and one around the $z$ axis. By Noether's theorem, if the Lagrangian is invariant under a one parameter differentiable transformation $g^\phi$, then there is a corresponding conserved quantity:
$$I(\mathbf{q}, \mathbf{\dot q}) = \frac{\partial L}{\partial \mathbf{\dot q}}\frac{dg^{\phi}(\mathbf{q})}{d\phi}|_{\phi=0}$$
You can find a proof of this theorem on pages 88-89 of Mathematical Methods of Classical Mechanics by V.I. Arnold. We have the good fortune of being in possession of three one parameter differentiable transformations under which $L$ is invariant, one for each axis of rotation:
$$g_x^{\phi} = \left(
\begin{array}{ccc}
1 & 0 & 0 \\
0 & \cos (\phi ) & -\sin (\phi ) \\
0 & \sin (\phi ) & \cos (\phi ) \\
\end{array}
\right) $$
$$g_y^{\phi} = \left(
\begin{array}{ccc}
\cos (\phi ) & 0 & \sin (\phi ) \\
0 & 1 & 0 \\
-\sin (\phi ) & 0 & \cos (\phi ) \\
\end{array}
\right)$$
$$g_z^{\phi} = \left(
\begin{array}{ccc}
\cos (\phi ) & -\sin (\phi ) & 0 \\
\sin (\phi ) & \cos (\phi ) & 0 \\
0 & 0 & 1 \\
\end{array}
\right)$$
Which gives us three more conserved quantities:
$$I_x(\mathbf{r}, \dot{\mathbf{r}}) = \dot z y - \dot y z $$
$$I_y(\mathbf{r}, \dot{\mathbf{r}}) = \dot x z - \dot z x$$
$$I_z(\mathbf{r}, \dot{\mathbf{r}}) = \dot y x - \dot x y$$
And these are the components of the angular momentum vector $\mathbf{L}$. Now we can rotate the axes so that $z$ points in the $\mathbf{L}$ direction, confining motion to the new $(x,y)$ plane. We actually lose two conserved quantities doing this, but we still have $\dot y x - \dot x y $ which I will just refer to as $L$ from here on out. Our now Hamiltonian now is just down to two degrees of freedom:
$$H = \frac{m_1m_2}{2(m_1+m_2)}(\dot x^2 +\dot y^2) + \frac{Gm_1m_2}{\sqrt{x^2+y^2}}\quad and \quad L = x\dot y- \dot x y$$
Or in polar coordinates,
$$H = \frac{m_1m_2}{2(m_1+m_2)}(\dot r^2 +r^2\dot \theta^2) + \frac{Gm_1m_2}{r}\quad and \quad L = r^2\dot\theta$$
And we can use $L$ to get rid of the second-to-last degree of freedom so we only have to deal with a first order ODE:
$$H = \frac{m_1m_2}{2(m_1+m_2)}(\dot r^2 +\frac{L^2}{r^2}) + \frac{Gm_1m_2}{r}$$
$$ \rightarrow\quad\frac{dr}{dt} = \sqrt{2(m_1+m_2)\left(\frac{H}{m_1m_2}-\frac{G}{r}\right)-\frac{L^2}{r^2}}$$
which can be integrated to get a relationship between $t$ and our last degree of freedom.
A few remarks. First, you may have noticed that the three angular momentum constants only let us reduce the number of degrees of freedom by two - this is because they are not functionally independent, which is the same as saying they do not form independent coordinates on the configuration space, which is the same as saying $\frac{\partial(I_x,I_y,I_z)}{\partial(x,y,z)} = 0$. This ends up not being a problem because with the Hamiltonian, we still have enough constants to reduce the problem to a first order ODE in a single degree of freedom.
The second thing is that there are actually other conserved quantities, namely the components of the eccentricity vector a.k.a the Laplace-Runge-Lenz vector which other commenters have mentioned. Unlike the constants of motion I've used above, the tactic of switching coordinate systems which reflect/reveal the symmetries of the system won't help you discover the LRL vector, at least not using simple coordinate transformations. That said, you can use the Hamilton-Jacobi method to solve this problem, in which case you'll find that the Hamilton-Jacobi PDE is solvable by separation of variables in more than one coordinate system (parabolic and spherical) which is how you can end up with more constants of motion than degrees of freedom. Evidently there is also a way to catch the LRL conserved quantities using a more involved treatment of Noether's theorem, but I haven't looked into it enough to comment on it. In any case, the LRL conserved vector unfortunately does not exist with different potentials or with 2+ bodies, whereas the conserved quantities outlined above are also conserved quantities of the $N$ body problem. Hope this helps.
