I would like to know how exactly the equations of motion in the Lorenz
gauge removes the second degree of freedom.
In the Lorenz 'gauge', we have
$$\Box A^{\mu} = \mu_0j^{\mu}$$
If $A^{\mu}$ is a solution, then so is $A^{\mu} + N\epsilon^{\mu}e^{-ik\cdot x}$ if
$$\Box (N\epsilon^{\mu}e^{-ik\cdot x}) = 0$$
Consistency with the Lorenz condition
$$\partial_{\mu}A^{\mu}=0$$
requires
$$k \cdot \epsilon = 0 $$
and consistency with the equation of motion requires that
$$k^2 = k \cdot k = 0$$
But this means that if a polarization four-vector $\epsilon$ satisfies the condition
$$k \cdot \epsilon = 0$$
then $\epsilon' = \epsilon + \alpha k$ also satisfies this condition
$$k \cdot \epsilon' = k \cdot (\epsilon + \alpha k) = k \cdot \epsilon + \alpha k^2 = 0$$
This means we can choose an $\epsilon^{\mu}$ such that $\epsilon^0 = 0$ and then the Lorenz condition implies that the wave and polarization 3-vectors are orthogonal
$$\vec k \cdot \vec\epsilon = 0$$
so there are only two independent polarization vectors (for freely propagating wave solutions).
To summarize, the Lorenz condition implies that the wave and polarization four-vectors are (Minkowski) orthogonal leaving three polarization degrees of freedom.
The equation of motion implies that the wave four-vector is null. Since null-vectors are self-orthogonal ($k^2 = 0$), we are left with two physical polarization degrees of freedom.
Comment to the question (v1): It seems OP is conflating, on one hand, a gauge transformation
$$ \tilde{A}_{\mu} ~=~ A_{\mu} +d_{\mu}\Lambda $$
with, on the other hand, a gauge-fixing condition, i.e. choosing a gauge, such e.g., Lorenz gauge, Coulomb gauge, axial gauge, temporal gauge, etc.
A gauge transformation can e.g. go between two gauge-fixing conditions. More generally, gauge transformations run along gauge orbits. Ideally a gauge-fixing condition intersects all gauge orbits exactly once.
Mathematically, depending on the topology of spacetime, it is often a non-trivial issue whether such a gauge-fixing condition is globally well-defined and uniquely specifies the gauge-field, cf. e.g. the Gribov problem. Existence and uniqueness of solutions to gauge-fixing conditions is the topic of several Phys.SE posts, see e.g. this and this Phys.SE posts.
Best Answer
The correct Gauge transformation formula should be $$\begin{aligned} \mathbf A &\mapsto \mathbf A + \nabla \lambda \\ \mathbf V &\mapsto V - \frac{\partial\lambda}{\partial t}, \end{aligned} $$ not something with "gradL/dt". The Coulomb gauge requires $\nabla\cdot\mathbf A=0$, not "rotA = 0". The Lorenz gauge requires $\nabla\cdot\mathbf A + \frac1{c^2}\frac{\partial V}{\partial t}=0$, not "gradA+1/c^2 dV/dt".
The Coulomb gauge can be chosen by solving the Poisson equation $$ \nabla^2 \lambda = -\nabla\cdot\mathbf A$$
The Lorenz gauge can be chosen by solving the inhomogeneous wave equation $$ \nabla^2 \lambda - \frac1{c^2}\frac{\partial^2\lambda}{\partial t^2} = -\nabla\cdot\mathbf{A} - \frac1{c^2}\frac{\partial{V}}{\partial{t}}$$
(Substitute the transformed potentials into the conditions to get the PDEs)
Existence of solutions of these PDEs are guaranteed as long as the source terms (stuff on the RHS) are "well-behaved" (e.g. $\nabla\cdot\mathbf A$ should grow slower than $1/r$ in the Poisson equation)