When studying the representation of the group $SO(3)$ on the space $L^2(\mathbf{S}^2)$, we have the spherical functions as its basis:
$$L^2(\mathbf{S}^2) = \text{span} \left \{ Y^l_m, l \in \mathbf{N}^+, -l \leqslant m \leqslant l \right \}$$ $$f(\theta,\varphi)=\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell f_\ell^m \, Y_\ell^m(\theta,\varphi)$$
But then suddenly the laplacian comes from nowhere and we says that its eigenfunctions are the spherical functions. I know it has important applications in physics, but from the group-theoretic point of view, is it just an operator like any other else?
Good question.
The Laplace-Beltrami on $\mathbb S^n$ is not "any" operator, of course. It is related to the group $SO(n)$, like the Euclidean Laplacian is related to the group of isometries of $\mathbb R^n$. The precise mathematical notion that captures this should be that of Casimir operator, but that is something that I do not understand, I must admit.
So I will give you my down-to-earth point of view. The starting point is Schur's lemma of representation theory, according to which, if $\phi\colon V\to V$ is an intertwining map of an irreducible representation of a group, then $\phi$ is a scalar multiple of the identity operator.
And that's precisely what the various Laplacians do: the Laplace-Beltrami commutes with $SO(n)$, therefore the spherical harmonics must be its eigenfunctions, because they are the basis functions of the irreps of $SO(n)$. Similarly, the Euclidean Laplacian commutes with translations, so the plane waves $e_\mathbf{p}(\mathbf{x}):=\exp(i\, \mathbf{x}\cdot\mathbf{ p})$ must be its eigenfunctions.
EDIT (after comments). I have deliberately been sloppy here, because there's a number of technical points that I do not want to address (mostly because I can't, to be honest).
I am sweeping something under the rug. Indeed, I am acting as if $\Delta$ were defined on $L^2(X)$, which is not the case: there exist functions whose Laplacian is not defined, or it is not an element of $L^2(X)$. This is especially bad in the case $X=\mathbb R^n$, in which case even the "eigenfunctions" that we discover are not eigenfunctions at all. I am talking about the plane waves $e_\mathbf{p}$, whose $L^2$ norm is $\infty$. There is a lot of theory to study to fully understand how to go around these technical difficulties.
But, whatever, I think that what we did in this post is enough as an illuminating guiding principle.