To prove explicitely that the Liénard-Wiechert potentials satisfy the Lorenz gauge, one has to find the time derivative and gradient of the retarded time. In this Wikipedia article it's calculated as follows:

Here $$\vec{r}$$ is the position vector of the observer, $$\vec{r}_s$$ the trayectory of the charged particle (evaluated at the retarded time) and $$\vec{n}=\frac{\vec{r}-\vec{r}_s|}{|\vec{r}-\vec{r}_s|}$$. Where I'm having some trouble is for the gradient, I don't understand where the $$\vec{n}$$ term comes from.

It took me a couple of steps of carefully taking the derivatives, but I can show it here. I note the source position with $$\textbf{r}'$$ instead of $$\textbf{r}_s$$ to avoid confusion with indices. Firstly the time derivative,
\begin{align} \frac{\partial}{\partial t} t &= \frac{\partial}{\partial t} \left(t_r + \frac{1}{c}|\textbf{r} - \textbf{r}'(t_r)|\right) \\ 1 &= \frac{\partial t_r}{\partial t} + \frac{1}{c}\frac{\partial}{\partial t}\sqrt{(\textbf{r} - \textbf{r}'(t_r))\cdot(\textbf{r} - \textbf{r}'(t_r))} \\ 1 &= \frac{\partial t_r}{\partial t} + \frac{1}{c}\frac{1}{2\sqrt{(\textbf{r} - \textbf{r}'(t_r))\cdot(\textbf{r} - \textbf{r}'(t_r))}} \frac{\partial}{\partial t}(\textbf{r} - \textbf{r}'(t_r))\cdot(\textbf{r} - \textbf{r}'(t_r))\\ 1 &= \frac{\partial t_r}{\partial t} + \frac{1}{c}\frac{\textbf{r} - \textbf{r}'(t_r)}{\sqrt{(\textbf{r} - \textbf{r}'(t_r))\cdot(\textbf{r} - \textbf{r}'(t_r))}} \cdot \left(\frac{\partial}{\partial t}\textbf{r} - \frac{\partial}{\partial t}\textbf{r}'(t_r) \right)\\ &\text{Note that }\frac{\partial}{\partial t}\textbf{r}=0,\\ 1 &= \frac{\partial t_r}{\partial t} - \frac{\hat{\textbf{n}}}{c}\cdot \frac{\partial}{\partial t}\textbf{r}'(t_r)\\ 1 &= \frac{\partial t_r}{\partial t} - \frac{\hat{\textbf{n}}}{c}\cdot \frac{\partial}{\partial t_r}\textbf{r}'(t_r)\frac{\partial t_r}{\partial t}\\ 1 &= \frac{\partial t_r}{\partial t} \left(1 - \frac{\hat{\textbf{n}}\cdot\textbf{v}'}{c} \right) \\ \frac{\partial t_r}{\partial t} &= \frac{1}{\left(1 - \frac{\hat{\textbf{n}}\cdot\textbf{v}'}{c} \right)}. \\ \end{align} Let's look at one component of the gradient, where repeated indices are summed over and $$\partial_i = \frac{\partial}{\partial x_i}$$, \begin{align} \partial_i t &= \partial_i t_r + \frac{1}{c}\partial_i |\textbf{r} - \textbf{r}'| \\ \partial_i t_r &= -\frac{1}{c} \partial_i \sqrt{(r_j - r'_{j}(t_r))(r_j - r'_j(t_r))}\\ \partial_i t_r &= -\frac{1}{c}\frac{(r_j - r'_{j}(t_r))}{|\textbf{r} - \textbf{r}'(t_r)|} \partial_i (r_j - r'_{j}(t_r))\\ &\text{Now we have the fact that in cartesian coordinates }\partial_i r_j = \delta_{ij}\\ \partial_i t_r &= -\frac{1}{c}\frac{(r_j - r'_{j}(t_r))}{|\textbf{r} - \textbf{r}'(t_r)|}(\delta_{ij} - \partial_i r'_{j}(t_r))\\ \partial_i t_r &= -\frac{1}{c}\frac{(r_i - r'_{i}(t_r))}{|\textbf{r} - \textbf{r}'(t_r)|} + \frac{1}{c}\frac{(r_j - r'_{j}(t_r))}{|\textbf{r} - \textbf{r}'(t_r)|}\frac{\partial r'_{j}(t_r)}{\partial t_r}\partial_i t_r\\ \partial_i t_r &= -\frac{1}{c}\frac{(r_i - r'_{i}(t_r))}{|\textbf{r} - \textbf{r}'(t_r)|} + \frac{1}{c}\frac{\textbf{r} - \textbf{r}'(t_r)}{|\textbf{r} - \textbf{r}'(t_r)|} \cdot \textbf{v}' \partial_i t_r\\ \end{align}
Now we can go back to dealing with the full vector, \begin{align} \nabla t_r &= -\frac{1}{c}\frac{\textbf{r} - \textbf{r}'(t_r)}{|\textbf{r} - \textbf{r}'(t_r)|} + \frac{1}{c} \left( \frac{\textbf{r} - \textbf{r}'(t_r)}{|\textbf{r} - \textbf{r}'(t_r)|} \cdot \textbf{v}' \right) \nabla t_r\\ \nabla t_r &= -\frac{1}{c}\hat{\textbf{n}} + \frac{1}{c} \left( \hat{\textbf{n}} \cdot \textbf{v}' \right) \nabla t_r\\ \nabla t_r &= -\frac{\hat{\textbf{n}}}{c\left(1 - \frac{\hat{\textbf{n}}\cdot\textbf{v}'}{c} \right)}. \\ \end{align}
As a side note, you can think of the spacetime derivative of $$t_r$$ as being a light-like covector meaning $$(\partial^\mu t_r)(\partial_\mu t_r) = 0$$. This gives, $$\frac{1}{c^2}\partial^0 t_r \partial_0 t_r = \nabla t_r \cdot \nabla t_r.$$ The fact that it is a covector describes the relative minus sign for the spatial component, and we can see that, $$\frac{1}{c}\partial_0 t_r = -|\nabla t_r|$$ This means that the only freedom left is a unit vector. The natural choice would be the normal to the hypersurface $$\hat{\textbf{n}}$$. Now you only need to work out the time derivative and the spatial derivative is fully determined.