Theorem. Let $n > 2$ be a positive integer. Assume that in the space
$\mathbb{R}^n$ (or possibly in $\mathbb{R}^n \setminus \{0\}$) we
have a vector field $g(r)$ with the property that for any $T \in
\text{SO}(n)$
$$T\,g(r) = g(Tr)$$
Then $g(r)$ is radial, i.e. there exists a scalar function $w :
\mathbb{R}_{+} \to \mathbb{R}$ such that
$$g(r) = w(|r|)\, r$$
Furthermore, if we define the function $$V(r) = W(|r|) \,\, \text{ where } \,\, W(\lambda) = \int_{\lambda_0}^{\lambda} \,s\, w(s) \, ds$$
then for any $r\in \mathbb{R}^n \setminus \{0\} $
$$\nabla \, V(r) = g(r)$$ i.e. all rotationally-invariant vector fields on $\mathbb{R}^n \setminus \{0\}$ are both radial and potential fields.
$$ $$
Proof: Let $r \in \mathbb{R}$ be an arbitrary non-zero vector. Define the satbilizer of $r$ from the rotation group $\text{SO}(n)$
$$\text{Stab}(r) = \{\, T \in \,\text{SO}(n) \, : \, T\,r = r\}$$ Furthermore, denote by $L$ the orthogonal complement of the vector $r$, which is by definition
$$L = \{\, v \in \mathbb{R}^n \, : \, v \cdot r = 0 \,\}$$
Then $\dim L = n-1$ and it is $\text{Stab}(r)$ invariant, i.e. for any $v \in \mathbb{R}^n$ such that $v \cdot r = 0$ and for any $T \in \text{Stab}(r)$
\begin{align*}
T\,v \cdot r = T\,v \cdot T\, r = v \cdot r = 0
\end{align*}
so $T\,v \in L$. Decompose the space $$\mathbb{R}^n = L \oplus \mathbb{R} \,r$$
Then $v \in \mathbb{R}^n$ decomposes uniquely as $v = v_L + \lambda \,r$, with $v_L \in L$.
For any $T \in \text{Stab}(r)$
$$T\, v = T(v_L + \lambda \,r) = T\, v_L + \lambda \,T \, r = w_L + \lambda \, r$$ where $w_L = T\,v_L \in L$.
Since $\text{Stab}(r)$ is a subgorup of $\text{SO}(n)$ and as such preserves the dot product on $\mathbb{R}^n$, when restricted to $L$ it preserves
the dot product on the $n-1$ dimensional subspace $L$. Moreover, if you take any linear transformation $T_L$ on $L$ that preserves the dot product,
it can be extended to $T \in \text{SO}(n)$ as $$T \,v = T\, (v_L + \lambda\, r) = T_L\, v_L + \lambda\, r$$ Thus, $\text{Stab}(r)$ is the full
rotation group of the $n-1$ dimensional subspace $L$ and is therefore isomorphic to $\text{SO}(n-1)$. Form the properties of $\text{SO}(n-1)$ we know
that if $v \in L$ is any non-zero vector, then its orbit under $\text{Stab}(r)$ is an $n-2$ dimensional sphere, so
there is always a transformation $T \in \text{Stab}(r)$ such that $T \,v \neq v$.
Now, let us focus on the $\text{SO}(n)$-invariant vector field $g(r)$. Take any non-zero vector $r \in \mathbb{R}^n$. Split the space
$$\mathbb{R}^n = L \oplus \mathbb{R} \,r$$ as explained above. Then $g(r)$ decomposes uniquely into
$$g(r) = g_L + w\, r$$
where $g_L \in L$ and $w \in \mathbb{R}$. For any $T \in \text{Stab}(r)$
$$T\,g(r) = g(T\,r) = g(r)$$ which in decomposed form translates into
$$T\,g(r) = T(g_L + w\, r) = T\,g_L + w\, T\,r = T\,g_L + w\, r = g_L + w\,r$$ and when we cancel out the term $w\,r$ from both sides of the last identity,
we find out that for any $T \in \text{Stab}(r)$
$$T\,g_L = g_L$$
But this is possible if and only if $g_L = 0$, because as I already mentioned before,
for a non-zero vector $g_L$ from $L$, there is always a transformation $T \in \text{Stab}(r)$ such that $T \,g_L \neq g_L$
(because the orbit is a proper $n-2$ dimensional sphere).
Hence, for any non-zero $r \in \mathbb{R}^n$ there exists a real number $u(r) \in \mathbb{R}$, that varies with respect to $r$,
such that $$g(r) = u(r)\, r$$ i.e. the vector field is radial.
Now, by invariance, for any $T \in \text{SO}(n)$ and for any non-zero $r \in \mathbb{R}^n$,
$$T (u(r)\, r) = u(r)\, T\,r = T\,g(r) = g(T\,r) = u(T\,r) \, T\,r$$ which, because of the identity $u(r)\, T\,r = = u(T\,r) \, T\,r$,
is possible if and only if $u(r) = u(T\,r)$.
Fix one unit vector $r_0 \in \mathbb{R}^n$ ($|r_0| = 1$). Take any other non-zero vector $r \in \mathbb{R}^n$. Both vectors $\frac{r}{|r|}$
and $r_0$ lie on the unit sphere and since $\text{SO}(n)$ acts transitively on the unit sphere (in fact with stabilizers always isomorphic to $\text{SO}(n-1)\,$)
there exists an orthogonal transformation $T \in \text{SO}(n)$ such that $$T \, \frac{r}{|r|} = r_0$$ Now do the following:
\begin{align}
u(r) =& u(T\,r) \\
u(r)\, =& \, u\left(\,T \left(|r| \, \frac{r}{|r|}\right)\,\right) \\
u(r) \, =&\, u\left(|r| \, T \, \frac{r}{|r|}\right) \\
u(r) \, =&\, u\big(\,|r| \, r_0\,\big)
\end{align}
It follows from the last identity that there exists a single-variable function $$w : \mathbb{R}_{+} \to \mathbb{R}$$
$$w(\lambda) = u(\lambda \, r_0)$$
such that
$$u(r) = w(|r|)$$
All of this allows us to conclude that for $n > 2$, any $\text{SO}(n)$-invariant vector field $g(r)$ on $\mathbb{R}^n$ is radial and can be expressed as
$$g(r) = w(|r|)\, r $$ for some suitable function $w : \mathbb{R}_{+} \to \mathbb{R}$. If we set
$$W(\lambda) = \int_{\lambda_0}^{\lambda} \, s\, w(s)\,ds$$
$$V(r) = W(|r|)$$ then $\frac{d}{d\lambda} W(\lambda) = W'(\lambda) = \lambda\, w(\lambda)$ and hence by the chain rule
$$\nabla\, V(r) = \nabla \, W(|r|) = W'(|r|) \, \frac{r}{|r|} = |r|\, w(|r|) \, \frac{r}{|r|} = w(|r|)\, r = g(r)$$ i.e.
the rotationally-invariant vector field $g(r)$ is potential.
$$ $$
$$ $$
Older Version.
Assume that for an odd dimensional space $\mathbb{R}^n$ we have a vector field $g(r)$ with the property that for any $T \in \text{SO}(n)$
$$T\,g(r) = g(Tr)$$
Theorem. Let $n$ be odd. For every non-zero vector $r \in \mathbb{R}^n$ there exists an orthogonal transformation $T \in \text{SO}(n)$ such that $Tr = r$ and the space of all vectors such that $Tv = v$ is one dimensional, spanned by $r$. In other words, for any vector $r \in \mathbb{R}^n$ there exists an orthogonal transformation $T \in \text{SO}(n)$ such that $r$ is an egienvector of $T$ with eigenvalue $1$, where the multiplicity of the eigenvalue is $1$.
Proof: Take an arbitrary non-zero vector $r \in \mathbb{R}^n$ and define the unit vector $\hat{r} = \frac{r}{|r|}$. Denote by $L$ the $n-1$ dimensional vector subspace of $\mathbb{R}^n$ orthogonal to $\hat{r}$. Since $L$ is non-trivial, choose a unit vector $\hat{r}_2 \in L$. Then $\hat{r}$ and $\hat{r}_2$ are pairwise orthogonal. Continue this procedure. Assume you have selected $k$ unit vectors $\hat{r}, \, \hat{r}_2, \, ..., \, \hat{r}_k$, that are pairwise orthogonal. Then, there exists a unique $n-k$ dimensional vector subspace $L_k$ which is the orthogonal complement of the span of $\hat{r}, \, \hat{r}_2, \, ..., \, \hat{r}_k$. If $k < n$ then $L_k$ is non-trivial and you can find aunit vector $\hat{r}_{k+1} \in L_k$. Thus, you have found $k+1$ unit vectors $\hat{r}, \, \hat{r}_2, \, ..., \, \hat{r}_k, \, \hat{r}_{k+1}$, that are pairwise orthogonal. Continue like this until $k = n$. You have constructed an orthonormal basis $\hat{r}, \, \hat{r}_2, \, ..., \, \hat{r}_n$ of $\mathbb{R}^n$. Think of this new vectors as column-vectors. Define the orthogonal matrix
$$U = \big[ \hat{r} \, \, \hat{r}_2 \,\, ... \, \hat{r}_n \big] \, \in \, \text{SO}(n)$$ If you denote by $e_1 = [1, \, 0, \, ..., \, 0]^T$ the first basis vector of $\mathbb{R}^n$, interpreted as a column-vector, then
$$\hat{r} = U \, e_1 \,\, \text{ and therefore } \,\, e_1 = U^T \, \hat{r}$$
If $n = 2\,m + 1$, construct the block-diagonal matrix
$$R = R(\theta_1,..., \theta_m) = \begin{bmatrix} 1 & & & &\\
& R_2(\theta_1) & & & &\\
& & R_2(\theta_2) & & \\& & & ... & \\ & & & & R_2(\theta_m)
\end{bmatrix}$$
where for each $j=1,..,m$
$$R_2(\theta_j) = \begin{bmatrix} \cos(\theta_j) & -\,\sin(\theta_j)\\
\sin(\theta_j) & \cos(\theta_j)
\end{bmatrix}$$ and $\theta_j \in (0, 2\pi)$. By construction, $$R(\theta_1,..., \theta_m)\, v = v \,\, \text{ if and only if } \,\, v = \lambda \,e_1$$, i.e. this matrix has eigenvalue $1$ with multiplicity $1$. The rest of the eigenvalues are true complex unitary numbers. Finally, define the orthogonal transformation
$$T = U \, R \, U^T \in \text{SO}(n)$$
By construction, $$T \, r = T \, |r| \, \frac{r}{|r|} = |r|\, T \,\hat{r} =
|r|\, U\, R \, U^T \,\hat{r} = |r|\, U\, R\, e_1 = |r|\, U\, e_1 = |r|\, \hat{r} = r$$ Moreover, since $T$ and $R$ are conjugate matrices, they have the same eigenvalues, including multiplicities. Therefore $T$ has eigenvalue $1$ of multiplicity $1$ and eigenvector $r$. End of proof
Now, you can apply this theorem to any $\text{SO}(n)$-invariant vector field on $\mathbb{R}^n \setminus \{0\}$, i.e. $$g : \mathbb{R}^n \setminus \{0\} \to \mathbb{R}^n$$
$$g(T\,r) = Tg(r) \,\,\text{ for any } \,\, T \in \text{SO}(n)$$ For an arbitrary non-zero vector $r \in \mathbb{R}^n$, pick one orthogonal transfromation $T_r \in \text{SO}(n)$ with the property described in the theorem above. Hence $T_r\, r = r$. Furthermore
$$T_rg(r) = g(T_rr) = g(r)$$ which means that $g(r)$ is also an eigenvector of $T_r$ with egienvalue $1$. Since the eigenspace of $T_r$ that corresponds to the eigenvaule $1$ is one dimensional and spanned by $r$, the vector $g(r)$ must be a multiple of $r$, i.e. there exists a scalar $\lambda(r) \in \mathbb{R}$ such that $$g(r) = \lambda(r) \, r$$ As you can see, the vector field $g(r) = \lambda(r) \, r$ is radial.
What happens if you take $n=2$? Can you find a counterexample? (I think it is fairly straightforward)
Now, there is another approach to your question, in particular that if you define the potential function
$$V(r) = G \, \int_{\Sigma} \, \frac{\rho(s)}{|r - s|} \, ds$$ and the function $\rho(s)$ is nice enough so that we can differentiate with respect to $r$ under the integral, then
$$g(r) = \nabla \,V(r) = - \, G \, \int_{\Sigma} \, \rho(s)\,\frac{(r-s)\,\,}{|r - s|^3} \, ds$$
So basically, your $\text{SO}(n)$-invariant vector field is in fact a potential field and comes from a potential function $V(r)$ which is also $\text{SO}(n)$-invariant, i.e. for any $T \in \text{SO}(n)$
$$V(Tr) = V(r)$$ Fix one unit vector $r_0 \in \mathbb{R}^n$ ($|r_0| = 1$). Now take any other non-zero vector $r \in \mathbb{R}^n$. Both vectors $\frac{r}{|r|}$
and $r_0$ lie on the unit sphere and since $\text{SO}(n)$ acts transitively on the unit sphere (in fact with stabilizers always isomorphic to $\text{SO}(n-1)\,$) there exists an orthogonal transformation $T \in \text{SO}(n)$ such that $$T \, \frac{r}{|r|} = r_0$$ Now do the following:
\begin{align}
V(r) =& V(T\,r) \\
V(r)\, =& \, V\left(\,T \left(|r| \, \frac{r}{|r|}\right)\,\right) \\
V(r) \, =&\, V\left(|r| \, T \, \frac{r}{|r|}\right) \\
V(r) \, =&\, V\big(\,|r| \, r_0\,\big)
\end{align}
It follows from the last identity that there exists a single-variable function $$W : \mathbb{R}_{+} \to \mathbb{R}$$
$$W(\lambda) = V(\lambda \, r_0)$$
such that
$$V(r) = W(|r|)$$ Then, by the chain rule, the gradient of $V$ with respect to $r$ is
$$\nabla \, V(r) = \nabla \, W(|r|) = \frac{dW}{d\lambda}(|r|)\, \frac{r}{|r|} = W'(|r|)\, \frac{r}{|r|}$$
i.e. the vector field $g(r) = \nabla \, V(r) = W'(|r|)\, \frac{r}{|r|}$ is radial.
