Dirac Delta and Sloppy Notation
The notation $$\delta_{x_0}(\phi) = \phi(x_0)$$ makes sense to me. It's clear that $\delta_{x_0}$ is a linear$^{[a]}$ functional, i.e. it takes a function $\phi$ as an input and produces a number $\phi(x_0)$ as output.
Suppose we'd like an integral representation of $\delta_{x_0}$. One way to do it is with a limit, e.g.
$$\delta_{x_0}(\phi) \equiv \phi(x_0) = \lim_{\sigma \rightarrow 0} \int_{-\infty}^\infty \underbrace{ \left[ \frac{1}{2 \pi \sigma^2} \exp \left( - \frac{(x-x_0)^2}{2 \sigma^2 }\right) \right]}_{\text{Approximate delta function}} \phi(x) \, dx \, .$$
The thing marked by the underbrace has integral 1, and as $\sigma \rightarrow 0$, it becomes thinner and taller. It's just annoying to write the limit and the Gaussian function in the integral, so physicists use a shorthand and write
$$\delta_{x_0}(\phi) \equiv \phi(x_0) = \int_{-\infty}^\infty \delta(x-x_0) \phi(x) \, dx \, . \tag{$\star$}$$
Here we're thinking of $\delta(x-x_0)$ as a function. Of course, as you say, it's not really function; there is no function that satisfies the integral equation $(\star)$. Indeed, $\delta(x-x_0)$ is better thought of as a limit of a function, as we wrote in the first integral equation.
Now let's answer the questions.
why does the $\delta(x−x_0)$ notation come in?
Well, we just showed that it more or less comes from the fact that, in a physicist's sense, $\delta$ stands for a limit of a real function, and that real function has argument $(x-x_0)$ (look at the Guassian in the first integral equation). That's a pretty weak answer though, but...
What mathematical power does one gain by representing the Dirac Delta as such?
...we can see why the $\delta(x-x_0)$ is powerful through an example. Consider the equation
$$ f(x_0) = \int_{-\infty}^\infty \delta(x-x_0) f(x) \, dx \, . \tag{$\star \star$}$$
Let's represent $f$ by its Fourier transform
$$f(x) = \int_{-\infty}^\infty \frac{dk}{2\pi} \tilde{f}(k) e^{ikx} \, . $$
Inserting the Fourier transform into $(\star \star)$ gives
\begin{align} f(x_0) &= \int_{-\infty}^\infty \int_{-\infty}^\infty dx \, \frac{dk}{2\pi} \delta(x-x_0) \tilde{f}(k) e^{ikx} \\ (\text{change variables }y\equiv x-x_0) \quad &= \int_{-\infty}^\infty \int_{-\infty}^\infty dy \, \frac{dk}{2\pi} \delta(y) \tilde{f}(k) e^{ik(y+x_0)} \\ &= \int_{-\infty}^\infty \frac{dk}{2\pi} \tilde{f}(k) e^{i k x_0} \\ &\equiv f(x_0) \, . \end{align}
The point here is that we can manipulate the argument of $\delta(x-x_0)$ e.g. in changes of variables just like any other function and it works. It's pretty obvious that this should work though, because as we pointed out above, you can always think of a $\delta$ as a limit of a sequence of ever-narrower real functions.
$[a]$: I say it's linear because
$$\delta_{x_0}(\phi + \theta) \equiv (\phi + \theta)(x_0) = \phi(x_0) + \theta(x_0) \equiv \delta_{x_0}(\phi) + \delta_{x_0}(\theta)\, .$$
\begin{equation} \beta(x_0) = \int_{-\infty}^{\infty}\alpha(x_0-x)\phi(x)dx\end{equation}
This is a convolution.
Convolutions are quite useful things to work with and show up all over the place.
If you know signal processing, the pointwise multiplication of a fourier transform of two functions is equal to the fourier transform of the convolution of two functions. And the fourier transform is its own inverse.
In signal processing, you use the fourier transform to extract the frequency components of an input function. Pointwise multiplication in the "frequency domain" corresponds to convolution in the time domain; so you can use convolution to do a "band pass" filter on a function (say, drop all components above a certain frequency).
This is all relatively academic, but the point is that convolutions are common and quite useful in a wide variety of places.
Now, if we imagine the Dirac Delta as a function we can convolve with another function such that when we do it, we extract the original function out.
\begin{equation} \int_{-\infty}^{\infty}\delta(x-x_0)\phi(x)dx = \phi(x_0)\end{equation}
You'll notice how this looks much like a convolution. Such forms are well studied and you can use other functions in place of the Dirac delta there and get reasonably interesting results.
For this to be "true", we just need delta to be a function that is zero everywhere except at 0, and its integral from negative epsilon to positive epsilon be 1. Then the math of convolution gives us that the result of the convolution is \phi(x_0).
\begin{equation} \int_{-\infty}^{\infty}\delta(x_0)\phi(x)dx = \beta(x_0)\end{equation}
Here, if we use the same imaginary Dirac delta, we end up with $\beta(x_0)$ equal to infinity at $x_0=0$, and equal to zero everywhere else.
This isn't $\phi(x_0)$ as we want.
The $x-x_0$ term is non-zero unless $x=x_0$, and zero if $x=x_0$. In a sense, this forces the entire "sampling" of $\phi$ to occur at $x=x_0$ and zero contribution to occur elsewhere.
To make this formal without using extended definitions of "function", we state that instead of the Dirac delta being one function, it is instead a limit of functions, and the convolution is the limit of the convolutions.
We take any ordered set of functions whose limit off 0 is 0, and whose integral from negative epsilon to positive epsilon limits to 1 from below for any epsilon.
For example, $\delta_i(x) = i/2$ if $|x|<1/i$ and 0 otherwise. Or a myriad of other options.
Now we take this:
\begin{equation} \int_{-\infty}^{\infty}\delta_i(x_0-x)\phi(x)dx = \beta_i(x_0)\end{equation}
This gives us a series of functions $\beta_i$, generated by a traditional convolution. The limit of the $\beta_i$ is $\phi$ regardless of what form $\delta_i$ take.1
Now, the equation you showed was $x-x_0$ instead of my convolution $x_0-x$. Well, the $\delta_i$ I used are symmetric around 0 -- so $\delta_i(x-x_0)=\delta_i(x_0-x)$.
In short, what is going on is that convolution is useful, the "pretend" Dirac delta treated as a convolution equation is often useful, and what is "really" going on can be interpreted as a limit of convolution of functions for which the equation stands in for.
Now this being mathematics, when we have a useful concept like Dirac delta, we go off and redefine what the symbols mean to make the Dirac delta a "generalized function" and make the notation we want to use "just work" without it "being shorthand for what is really going on", because all practical math is just shorthand anyhow.
Finally, $\delta(x-x_0)(\phi)$ just looks like sloppy notation. They are using $\delta$ as a kind of macro to me. I don't see the use of it, really.
I'd want $\delta_\phi(x_0)$ or $\delta(\phi)(x_0)$ or even $\delta_{x_0}(\phi)$ or $\delta(x_0)(\phi)$ or or $\delta * \phi$ equal to $\int_{-\infty}^{\infty}\delta_i(x_0-x)\phi(x)dx$ (where $*$ is convolution)
1 My fourier analysis is a bit rusty, I'd believe there are pathological functions that don't behave well when convoluted yet satisfy the requirements I placed on the $\delta_i$. So, I'll hedge my bets and say "sufficiently nice" functions instead.
As DanielSank pointed out, the notation $\delta(x-x_0)(\phi)$ is non-standard.
As to where "physicists' notation" originates, it comes from the representation of the delta distribution as a limit
$$ \delta_{x_0}[\phi] = \lim_{\epsilon\to0} \int \delta_\epsilon(x-x_0)\phi(x)\mathrm dx $$
of functions $\delta_\epsilon$ such as
$$ \delta_\epsilon(x) = \frac1\pi \frac\epsilon{x^2 + \epsilon^2} $$