For me, the gradient of a scalar field (say, in three dimensions) is simply (formally)
$\nabla f = \left(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y},\frac{\partial f}{\partial z} \right)$.
In which way do we need a metric?
But some people did tell me only on a Riemann manifold (with a metric) can we define the gradient.
On any manifold we can define the differential $df$ of a scalar $f$. The differential is a 1-form: something that eats vectors and spits out scalars, or even less formally, something with one down index. We have the following formula for the differential, $$df = \frac{\partial f}{\partial x^i} dx^i$$ (sum over $i$ implied). You can write it in index notation as $$(df)_\mu = \frac{\partial f}{\partial x^\mu}.$$
If by gradient you mean a vector field, then to make 1-forms into vector fields, you need something to "raise the index". That's where the metric comes in.
The main point of the gradient is to have the following equation: $$ df(v) = \langle \nabla f , v \rangle $$ for every vector $v$ of the tangent space where $\langle \cdot , \cdot \rangle$ is the metric.
In local coordinates $(x^{i})$, using the Einstein notation, the dot product of $v = v^{i}\frac{\partial }{\partial x^{i}}$ and $w = w^{j}\frac{\partial }{\partial x^{j}}$ can be written as: $$ \langle v , w \rangle = g_{ij}v^{i}w^{j} $$
The diferential of $f$ is $df = \frac{\partial f}{\partial x^{i}}dx^{i}$ (notice that is independent of the metric) so if we write $v = v^{j}\frac{\partial }{\partial x^{j}}$ we have: $$ df(v) = \frac{\partial f}{\partial x^{i}}dx^{i}(v) = \frac{\partial f}{\partial x^{i}}dx^{i}(v^{j}\frac{\partial }{\partial x^{j}}) = \frac{\partial f}{\partial x^{i}} v^{j} dx^{i}(\frac{\partial }{\partial x^{j}}) = \frac{\partial f}{\partial x^{i}} v^{j} \delta_{j}^{i} = \frac{\partial f}{\partial x^{i}} v^{i} $$
In the other hand, $ \langle \nabla f , v \rangle = g_{ij} (\nabla f)^{i}v^{j} $, where since we want $\nabla f$ being a vector we write $ \nabla f = (\nabla f)^{i}\frac{\partial }{\partial x^{i}}$.
We remark now that if we set $ (\nabla f)^{i} = g^{ij}\partial_{j}f$ (where $g^{ij}$ is the inverse of the metric and $\partial_{j}f = \frac{\partial f}{\partial x^{j}}$) we have: $$ \langle \nabla f , v \rangle = g_{ij} (\nabla f)^{i}v^{j} = g_{ij} g^{ik}\partial_{k}fv^{j} = \delta_{j}^{k}\partial_{k}fv^{j} = \partial_{j}fv^{j} = \frac{\partial f}{\partial x^{i}} v^{i} = df(v) $$
which is what we wanted.
Now it is clear that the gradient of $f$ is dependent of the metric since we use the inverse of the metric to define it.
It's all about making things coordinate-independent. In the orthonormal standard basis of Euclidean space $\{\mathbf e_i\}$, the gradient reads $$ \nabla f = \sum_i \frac{\partial f}{\partial x^i} \mathbf e_i $$ Now, choose a different set of coordinates $x'=\phi(x)$. Then, we have $$ f = f'\circ\phi $$ and thus $$ \frac{\partial f}{\partial x^i} = \sum_j \frac{\partial f'}{\partial x'^j} \frac{\partial\phi^j}{\partial x^i} $$ as well as $$ \mathbf e_i = \sum_k \frac{\partial\phi^k}{\partial x^i} \mathbf e'_k $$ which is somewhat non-obvious.
Putting these together, we end up with $$ \nabla f = \sum_{ijk} \frac{\partial\phi^j}{\partial x^i} \frac{\partial\phi^k}{\partial x^i} \frac{\partial f'}{\partial x'^j} \mathbf e'_k $$
This can be written in a form-invariant way $$ \nabla f = \sum_{ij} g^{ij} \frac{\partial f}{\partial x^i} \mathbf e_j = \sum_{ij} g'^{ij} \frac{\partial f'}{\partial x'^i} \mathbf e'_j $$ where $$ g^{ij} = \delta^{ij} $$ is the Euclidean metric on the cotangent space and $$ g'^{ij} = \sum_{kl} \frac{\partial\phi^i}{\partial x^k} \frac{\partial\phi^j}{\partial x^l} g^{kl} $$ its transform.
Note that in contrast to the gradient, the differential already is form-invariant $$ \mathrm df = \sum_i \frac{\partial f}{\partial x^i} \mathrm dx^i = \sum_i \frac{\partial f'}{\partial x'^i} \mathrm dx'^i $$ independent of any metric.
This leads to the coordinate-free definition of the gradient in terms of the metric $g$ on the tangent space $$ \nabla f\rfloor g = \mathrm df $$
$\mathbb{R}^3$ is a Riemannian manifold with metric $\text{diag}(1, 1, 1)$, and on this manifold with the specified metric the gradient takes the form you've given. With a more general manifold the metric-compatible derivative has extra terms given by the Christoffel symbols, but the gradient of a scalar field will still be given by the partial derivatives.
I find it simplest to define the total derivative first, since it is independent of the metric and noting that it is a linear operator mapping differentials to differentials: $$u = \phi(\mathbf{x}), \quad du = \phi'(\mathbf{x})[d\mathbf{x}]$$ Noting that the differential $d\mathbf{x}$ is simply a vector (difference in position vectors). The derivative, $\phi'$ at $\mathbf{x}$ is a linear functional, a mapping from a vector to a scalar. It thus is an element of the dual space $X^*$, and it is called a dual vector.
To define it in terms of components we invoke partial derivatives: $$\phi'(x): d\mathbf{x} \mapsto \frac{\partial \phi}{\partial x} dx + \frac{\partial \phi}{\partial y} dy + \ldots$$ In words: $\phi'$ is the linear functional mapping the basis vectors: $\hat{\imath}$ to $\partial_x \phi$ and $\hat{\jmath}$ to $\partial_y\phi$ and $\hat{k}$ to $\partial_z \phi$. (or something similar in other coordinates.)
While this derivative as an operator is not itself a vector but a dual vector, the Reisz Representation Theorem says that (in finite dimensions) any dual vector can be expressed by taking the inner (dot) product (which is the metric) with respect to a unique vector. $$ \xi \in X^*, \exists \mathbf{u}: \xi(\mathbf{v}) =\mathbf{u}\bullet \mathbf{v}$$
In particular the vector which does this for the total derivative functional is the gradient vector. $$ \phi'(\mathbf{x})[d\mathbf{x}] = \nabla\phi \bullet d\mathbf{x}$$
So the answer to your question is that to get from the (metric independent) derivative to the gradient we must invoke the metric. In component form (summing over repeated indices):
$$ \nabla \phi = g^{\mu\nu}\frac{\partial \phi}{\partial x^\mu} \mathbf{e}_\nu$$ The coordinates have raised indices to contract with the lower indices of the basis to which they are coefficients. Thus the partial derivatives have lowered indices which must be raised (by the metric) to again contract with basis vectors.
And finally, when you look at what happens under general change of bases and change of coordinate systems, the partial derivatives do not transform like components of a vector but the partial derivatives when contracted with the metric do so transform.
I'm not a physicist and not great with tensor technicalities but I did have a desire to understand this issue, having encountered the Fisher Information metric in the field of information geometry.
So here's a very informal perspective that connects this to a super basic linear algebra procedure. I think it's roughly correct in spirit and definitely helped me understand this.
This borrows some ideas from MathTheBeautiful's tensor calc series. I'd recommend his video on the covariant basis if you're not familiar with it.
Let $f$ be a function defined in some "raw geometric" coordinate-independent space. Let $\vec{R}$ denote a vector in this space. So $f(\vec{R})$ is the value of $f$ at the point in the space $\vec{R}$. We could think of the "true gradient" (maybe abusing notation) as $\frac{df}{d\vec{R}}$. This gradient is the vector such that the inner product of another vector $v$ in the space with $\frac{df}{d\vec{R}}$ gives you the directional derivative (i.e., the linear part of the increase of $f$ in the direction of $v$ with "speed" proportional to the magnitude of $v$).
Say we have some random coordinate system $\left[s, t, u \right]$. The covariant basis is the set of derivatives of a position vector $\vec{R}$ in the geometric space with respect to each coordinate.
So, at some point of evaluation in our space, the covariant basis will be $\frac{d\vec{R}}{ds}, \frac{d\vec{R}}{dt}, \frac{d\vec{R}}{du} $
The goal is to take the partials of our function $f$ w.r.t. our coordinates ($\frac{\partial f}{\partial s}, \frac{\partial f}{\partial t}, \frac{\partial f}{\partial u} $) that we started with, and use those to get coefficients for our (covariant) basis so that we can express the gradient in the basis.
By definition of the gradient (or you could kind of think about the chain rule), each partial derivative is the inner product of the corresponding basis vector with the gradient, $\frac{\partial f}{\partial s} = \frac{df}{d\vec{R}} \cdot \frac{d\vec{R}}{ds}$.
This makes sense. The partial w.r.t. $s$ is the directional derivative of $f$ in the direction of $\frac{d\vec{R}}{ds}$, since that's the direction you move when you increase $s$.
Now for what's probably really abusive notation. Let $B$ be the matrix with the covariant basis as column vectors. Then we can express the "vector" of partials as follows
$$ \left(\begin{matrix} \frac{\partial f}{\partial s} \\ \frac{\partial f}{\partial t} \\ \frac{\partial f}{\partial u} \end{matrix}\right) = B^T \frac{df}{d\vec{R}} $$
We're just matrix-multiplying $B^T$ by the gradient to get the "list" of inner products discussed above.
Let's suppose that we can express the gradient in our basis (our end goal after all). We don't actually know what the coefficients would be but we'll write them down as unknowns and try to recover the coefficients. Suppose
$$ \frac{df}{d\vec{R}} = a \frac{d\vec{R}}{ds} + b \frac{d\vec{R}}{dt} + c \frac{d\vec{R}}{du} = B \left(\begin{matrix} a \\ b \\ c \end{matrix}\right) $$
Now let's substitue the far-right expression into the above equation for the partials.
$$ \left(\begin{matrix} \frac{\partial f}{\partial s} \\ \frac{\partial f}{\partial t} \\ \frac{\partial f}{\partial u} \end{matrix}\right) = B^TB \left(\begin{matrix} a \\ b \\ c \end{matrix}\right) $$
Note $B^TB$ is the metric tensor. It's the matrix of all pairwise inner products of our covariant basis vectors. And it should be invertible. So to recover $a$, $b$, and $c$, we just multiply both sides by $(B^TB)^{-1}$.
$$ (B^TB)^{-1} \left(\begin{matrix} \frac{\partial f}{\partial s} \\ \frac{\partial f}{\partial t} \\ \frac{\partial f}{\partial u} \end{matrix}\right) = \left(\begin{matrix} a \\ b \\ c \end{matrix}\right) $$
The true gradient is some particular linear combination of our covariant basis vectors. The partial derivatives w.r.t. coordinates are each of the respective basis vectors dotted with this linear combination. That's where the metric tensor, the matrix of pairwise inner products of the covariant basis, $B^TB$ comes from. To recover the coefficients of this linear combination, we have to "invert off" the $B^TB$.
If by the word "gradient" you mean the 1-form whose components are $$ \partial_a \phi $$ then you do not need a metric to define a gradient. If by the word "gradient" you mean the associated vector field whose components are $$ g^{a \mu} \partial_\mu \phi $$ then you need a metric (or some other tool to map from cotangent space to tangent space) to define a gradient. Similar statements apply to $\nabla_a \phi$ and $\nabla^a \phi$.
(This repeats the answer by Robin Ekman; my aim is to clarify).
