Expanding on my comment: I will show that even looking only locally, this is impossible for dimension $n > 1$.

Define $\pi: TM \rightarrow T^*M$ as the "lowering the index" operator, taking a vector field $X \in TM$ to the covector field $\alpha \in T^*M$, where $\alpha(Y) = \langle X, Y\rangle$ for any vector field $Y$. This operator $\pi$ induces a natural linear operator from $\Gamma(TM)$ to $\Gamma (T^{*} (M)$. We again use $\pi$ to denote this induced operator. For the sake of notation, define $\{f, g\} = \langle \nabla f, \nabla g \rangle = (\nabla f)(g) = (\nabla g)(f)$

By definition, $\nabla f = \pi^{-1}(df)$. Therefore, for any gradient, we will have that $d(\pi(\nabla f)) = d(df) = 0$. We therefore only need to show that there are $f, g$ with $d(\pi([\nabla f, \nabla g])) \neq 0$. Note that this is in fact an if-and-only-if if we work locally.

Because the metric is nondegenerate, we can check what happens for vector fields that are gradients. We'll contract with one gradient at a time.

$\iota_{\nabla h}(d(\pi([\nabla f, \nabla g]))) = \mathcal{L}_{\nabla h}(\pi([\nabla f, \nabla g])) - d(\pi([\nabla f, \nabla g])(\nabla h))$

The second term is equal to $d(\langle [\nabla f, \nabla g], \nabla h \rangle)$. By the definition of $\nabla h$, this is equal to $d(\{f, \{g, h\}\} - \{g, \{f, h\}\})$.

$\iota_{\nabla j} \iota_{\nabla h}(d(\pi([\nabla f,\nabla g]))) = \iota_{\nabla j}(\mathcal{L}_{\nabla h}(\pi([\nabla f, \nabla g])) - d(\{f, \{g, h\}\} - \{g, \{f, h\}\}))$

$= \mathcal{L}_{\nabla h}(\pi([\nabla f, \nabla g])(\nabla j)) - \iota_{\mathcal{L}_{\nabla h}(\nabla j)}(\pi([\nabla f, \nabla g])) - \{j, \{f, \{g, h\}\}\} + \{j, \{g, \{f, h\}\}\}$

$= - \langle [\nabla f, \nabla g], [\nabla h, \nabla j]\rangle + \{h, \{f, \{g, j\}\}\} - \{h, \{g, \{f, j\}\}\} - \{j, \{f, \{g, h\}\}\} + \{j, \{g, \{f, h\}\}\}$

So we want to check if this is true for every $f, g, h, j$:

$\langle [\nabla f, \nabla g], [\nabla h, \nabla j]\rangle = \{h, \{f, \{g, j\}\}\} - \{h, \{g, \{f, j\}\}\} - \{j, \{f, \{g, h\}\}\} + \{j, \{g, \{f, h\}\}\}$

Call the left side $A(f, g, h, j)$ and the right side $B(f, g, h, j)$. Assume that $A(f, g, h, j) = B(f, g, h, j)$ for all $f, g, h, j$. Then $A(f, g g', h, j) - g A(f, g', h, j) - g' A(f, g, h, j) = B(f, g g', h, j) - g B(f, g', h, j) - g' B(f, g, h, j)$.

A quick calculation shows that $A(f, g g', h, j) - g A(f, g', h, j) - g' A(f, g, h, j) = \{f, g\} \{h, \{j, g'\}\} - \{f, g\} \{j, \{h, g'\}\} + \{f, g'\} \{h, \{j, g\}\} - \{f, g'\} \{j, \{h, g\}\}$.

On the other hand, $B(f, g g', h, j) - g B(f, g', h, j) - g' B(f, g, h, j) = \{f, g\} \{h, \{g', j\}\} + \{h, g\} \{f, \{g', j\}\} + \{g, j\} \{h, \{f, g'\}\} + \{f, g'\} \{h, \{g, j\}\} + \{h, g'\} \{f, \{g, j\}\} + \{g', j\} \{h, \{f, g\}\} - \{h, g\} \{g', \{f, j\}\} - \{h, g'\} \{g, \{f, j\}\} - (\{f, g\} \{j, \{g', h\}\} + \{j, g\} \{f, \{g', h\}\} + \{g, h\} \{j, \{f, g'\}\} + \{f, g'\} \{j, \{g, h\}\} + \{j, g'\} \{f, \{g, h\}\} + \{g', h\} \{j, \{f, g\}\} - \{j, g\} \{g', \{f, h\}\} - \{j, g'\} \{g, \{f, h\}\})$

So we want to check whether:

$\{h, g\}(- \{g', \{f, j\}\} - \{j, \{f, g'\}\} + \{f, \{g', j\}\}) + \{j, g\}(\{g', \{f, h\}\} + \{h, \{f, g'\}\} - \{f, \{g', h\}\}) + \{h, g'\}(- \{g, \{f, j\}\} - \{j, \{f, g\}\} + \{f, \{g, j\}\}) + \{j, g'\}(\{g, \{f, h\}\} + \{h, \{f, g\}\} - \{f, \{g, h\}\}) = 0$.

Using a similar trick by splitting $f$ to $f f'$, we get:

$2 \{h, g\}\{g', f\}\{f', j\} + \text{similar terms} = 0$, where "similar terms" refers to all terms gotten by switching g and g', switching f and f', and switching h and j and negating.

We've finally made it to a pointwise condition - which means that we only need to check this on a vector space with the standard metric. And it is false; take $h = g = g' = f = x, f' = j = y$. More generally, if we let $h, g, g', f, f', j = a_i x + b_i y$, the equation is:

$(a_1 a_2 + b_1 b_2)(a_3 a_4 + b_3 b_4)(a_5 a_6 + b_5 b_6) + \text{similar terms} = 0$.

The terms consisting of just $a$ or just $b$ cancel out; the only terms left either have 4 $a$s and 2 $b$s or vice versa. It's not hard to see that these terms can't cancel - the two $b$s (or the two $a$s) must remain together, which gives a unique product it can come from.

To summarize: We took the question, extracted a differential condition from it, extracted a pointwise condition from the differential condition, and showed that the pointwise condition could not be satisfied in dimension > 1.

ETA: To give a more explicit answer: choose coordinates $x, y$. Then:

$d(\pi([\nabla x^2, \nabla xy])) = d(\pi([2x \nabla x, x \nabla y + y \nabla x])) = d(\pi(2x \nabla y - 2y \nabla x + 2x^2 [\nabla x, \nabla y]))$

$= d(2x dy - 2y dx + 2x^2 [\nabla x, \nabla y]) = 4 dx \wedge dy + 4x dx \wedge \pi([\nabla x, \nabla y]) + 2x^2 d(\pi([\nabla x, \nabla y]))$

Just by checking at 0, we can see that this is nonzero - so $d(\pi([\nabla x^2, \nabla xy]))$ must be nonzero - so as long as $dx$ and $dy$ are linearly independent at some point where $x = 0$ (hence the need for $n > 1$), then at that point, $[\nabla x^2, \nabla xy]$ is not a gradient.