The ostensibly sad reality is that any such $\psi_{x_o}$ isn't really a position eigenfunction. Try acting on it with the position operator. You simply don't get $x_o$ times the wavefunction.
That said, there are important and meaningful ways in which it is nearly a position eigenfunction. Notice that such a Gaussian is extremely narrow for very small values of $a$. Thus you can approximate multiplication by $x$ (which is the action of the position operator in the position basis) by multiplication by $x_o$. In the only places it matters (where the wavefunction isn't exceedingly small), the position is very nearly $x_o$, so it doesn't do too much damage to just pretend like it really is $x_o$.
Indeed, when we take the limit $a \to 0$ it is an exact result. The trouble is that the state becomes unnormalizable in this case. Sure, it's a delta function so it's area is still just one, but the area under its absolute square--what we're really interested in--is infinite. This shouldn't bother you too much though, because these true position eigenstates (delta functions in the position basis) form a complete basis for the space of physical wavefucntions. Any such state $\psi(x)$ can be represented as a superposition of position eigenstates:
$$ \psi(x) = \int dx'\, \psi(x') \delta(x-x') = \int dx'\, \psi(x')\psi_{x}(x')$$
Where now $\psi_x(x')$ represents the true (unnormalizable) position eigenstate at position x.