I've learned a few things over the last couple of years and would like to share. Much of this has been said above, but I'd like to collect things in one place. Here are a few big reasons the Laplacian is important to Riemannian geometry:
$(1)$ Hodge Theory.
The Hodge Theorem states that on a compact oriented Riemannian manifold $(M,g)$, one has:
$$\{\alpha \in \Omega^k(M) \colon \Delta \alpha = 0\} \cong H^k_{\text{dR}}(M).$$
This is remarkable: the left-hand side is the space of harmonic $k$-forms, which a priori depends on the metric. But the right-hand side depends only on the homotopy type of the manifold! This provides a link between analysis and geometry on the left, and topology on the right.
This theorem has been generalized in several directions. For example: If $(M,g)$ is a compact Hermitian manifold, then
$$\{\alpha \in \Omega^{p,q}(M) \colon \Delta_{\overline{\partial}}\alpha = 0\} \cong H^{p,q}(M).$$
Again, the left side a priori depends on both the complex structure and the metric, whereas the right side depends only the complex structure.
As mentioned in original post, the Hodge Theorem can be coupled with Weitzenbock formulas to prove vanishing theorems. The classic example: If $\omega$ is a harmonic $1$-form and $X = \omega^\#$ is its dual vector field, then
$$\textstyle \Delta \frac{1}{2}|X|^2 = |\nabla X|^2 + \text{Ric}(X,X).$$
The Bochner method then shows that if $(M^n,g)$ is compact, oriented, and has $\text{Ric} \geq 0$, then the harmonic $1$-forms coincide with the parallel $1$-forms, and hence (by the Hodge Theorem) $b^1(M) \leq n$. If we had the stronger $\text{Ric} > 0$, then $b^1(M) = 0$. This is one example of many.
$(2)$ Spectral Geometry.
Given a linear operator, we should examine its eigenvalues (and eigenspaces).
Example: On an isometrically immersed surface in $\mathbb{R}^3$, the eigenvalues of the shape operator are the principal curvatures of the surface.
Example: On an oriented Riemannian $4$-manifold, the eigenspaces of the Hodge $\ast$ operator give the decomposition of $2$-forms $\Lambda^2 = \Lambda^2_+ \oplus \Lambda^2_-$ into self-dual and anti-self-dual parts.
The eigenvalues of the Laplacian provide invariants of the Riemannian manifold, and so encode geometric information. Now, I still don't have an intuitive grasp of this information (at the time of writing), but I've learned a few cool things.
First, on a compact Riemannian manifold $(M^n,g)$, the lowest eigenvalue $\lambda_1$ of $\Delta$ is related to the Ricci curvature. That is, if $\text{Ric} \geq k > 0$, then (Lichnerowicz)
$$\lambda_1 \geq \frac{n}{n-1}k$$
and this estimate is sharp (Obata): equality holds iff $(M^n,g)$ is the round sphere.
Second, again on a compact Riemannian manifold $(M^n,g)$, the sequence of eigenvalues grow according to an asymptotic formula (Weyl's Law) that depends only on the dimension $n$ and the volume of $(M,g)$.
In general, the eigenvalues cannot be computed explicitly (except in very special cases), and it is estimates on the eigenvalues that provide links to geometry. That said, the sequence of eigenvalues alone is not enough to completely determine the geometry (Milnor tori).
$(3)$ Euler-Lagrange Equations.
A significant part of Riemannian geometry is tied up with the calculus of variations -- that is, with various functionals and their critical points and extrema. The basic examples are the volume functional and energy functionals -- forming the theory of geodesics, minimal submanifolds, and harmonic maps -- but there are many others beyond these.
The critical points of a functional are described by "Euler-Lagrange equations," and these equations frequently involve the Laplacian. For example, if $F \colon (M,g) \to (N,h)$ is an isometric immersion, then the condition that $F(M) \subset N$ have zero mean curvature is the PDE system given in local coordinates by
$$\Delta_M F^j - \sum_{i,k} \sum_{\alpha, \beta} g^{\alpha \beta}(x) \Gamma^j_{ik}(F(x)) \frac{\partial F^i}{\partial x^\alpha} \frac{\partial F^k}{\partial x^\beta} = 0,$$
for $1 \leq j \leq \dim(N)$.
Remarks
(a) The Laplacian can be placed in the more general context of (squares of) Dirac operators. I'd like to say more about this some other time.
(b) The condition $\Delta \alpha = 0$ is equivalent to having both $d\alpha = 0$ and $d^*\alpha = 0$. This recasts the 2nd-order equation $\Delta \alpha = 0$ as a 1st-order system.
This "closed and co-closed" pairing seems to come up a lot. Examples include the Yang-Mills equation (for connections) and the integrability condition for $\text{G}_2$-structures.
(c) Even in the setting of flat $\mathbb{R}^n$, harmonic functions (solutions to $\Delta u = 0$) are pretty interesting. For one thing, solving $\Delta u = 0$ implies that the function $u$ enjoys many properties similar to holomorphic functions:
- Mean Value Property
- An integral representation formula
- Liouville-type theorems
- Maximum Principle
- Regularity results
This suggests that many of the nice properties of holomorphic functions in $\mathbb{C}$ are really artifacts of harmonicity.
For another thing, solutions of $\Delta u = 0$ are exactly the critical points of the Dirichlet energy functional. This links harmonic functions to the calculus of variations. More geometrically, it suggests a link between harmonic functions and minimal submanifolds. There are several results which realize this link, but I'll leave it at that.