Let $x,y,z$ be a system of coordinates with unit vectors respectively $\mathbf{u}_x$, $\mathbf{u}_y$, $\mathbf{u}_z$. Moreover, we have
$\nabla = \displaystyle \frac{\partial}{\partial x} \mathbf{u}_x + \frac{\partial}{\partial y} \mathbf{u}_y + \frac{\partial}{\partial z} \mathbf{u}_z$
Let $x',y',z'$ be a new system of coordinates: $x' = x$, $y' = y$, $z' = -z$.
But what about the unit vector and the Nabla operator?
I suppose the new unit vector are $\mathbf{u}_{x'} = \mathbf{u}_x$, $\mathbf{u}_{y'} = \mathbf{u}_x$, $\mathbf{u}_{z'} = -\mathbf{u}_z$.
Question 1 Is it always so? May we change the coordinates without changing the unit vectors?
We should then have
$\nabla' = \displaystyle \frac{\partial}{\partial x'} \mathbf{u}_{x'} + \frac{\partial}{\partial y'} \mathbf{u}_{y'} + \frac{\partial}{\partial z'} \mathbf{u}_{z'}$
But if I try to use the old system of coordinates I would obtain
$\nabla' = \displaystyle \frac{\partial}{\partial x} \mathbf{u}_{x} + \frac{\partial}{\partial y} \mathbf{u}_{y} + \frac{\partial}{\partial (-z)} (- \mathbf{u}_{z}) = \nabla$
Question 2 Is it correct? Or should I have expected another result?
First, you should get it out of your head that $\nabla$ is defined in terms of partial derivatives and unit vectors. This is not generally the case. Typically, one has
$$\nabla = g^1 \partial_1 + g^2 \partial_2 + g^3 \partial_3$$
where the vectors $g^1, g^2, g^3$ are not necessarily unit (or orthogonal, for that matter).
Now then, let
$$f(r) = f(x e_1 + y e_2 + z e_3) = r' = x' e_1 + y' e_2 + z' e_3$$
You can compute the Jacobian $\underline f$ in the usual manner:
$$\underline f(a) = (a \cdot \nabla) f \implies \underline f(e_1) = e_1, \underline f(e_2) = e_2, \underline f(e_3) = -e_3$$
From here you could compute the transpose $\overline f$, but this transpose is trivially equal to $\underline f$, so I won't bother going through the calculation process.
Now, using the chain rule, you can show that, if you have a scalar field $\phi(r)$ and another corresponding scalar field $\phi'(r') = (\phi' \circ f)(r)$, the chain rule tells us that
$$a \cdot \nabla \phi' = [a \cdot \nabla f] \cdot \nabla' \phi' = \underline f(a) \cdot \nabla' \phi' = a \cdot \overline f(\nabla') \phi'$$
Or, more succinctly,
$$\nabla = \overline f(\nabla')$$
Explicitly, what I've written is the idea that
$$e^1 \frac{\partial}{\partial x} + e^2 \frac{\partial}{\partial y} + e^3 \frac{\partial}{\partial z} = \overline f \left[ e^1 \frac{\partial}{\partial x'} + e^2 \frac{\partial}{\partial y'} + e^3 \frac{\partial}{\partial z'} \right] = e^1 \frac{\partial}{\partial x'} + e^2 \frac{\partial}{\partial y'} - e^3 \frac{\partial}{\partial z'}$$
This process is valid for arbitrary coordinate transformations; in particular, if you were to transform to spherical coordinates, it would correctly spit out non-unit vectors for the angular coordinates.
Secondly, I would attribute no geometric significance to the object $\nabla'$. It is merely a handy construct. In almost all cases, what you're really interested in is finding an expression for $\nabla$ (which never changes) merely in some different coordinate system.