You can also do this "topologically".  The idea is to separately calculate the (l-adic etale) Euler characteristics of the stack $Y_0(p)$ in $\mathrm{char}$ $0$ and in  $\mathrm{char}$ $p$, then see how they must relate.  Here we go:

Over $\mathbb{C}$, and hence over $\bar{\mathbb{Q}}_p$ as well, the Euler characteristic of $Y_0(p)$ is $(p+1)\times(-1/12)$, since $Y_0(p)$ is a $(p+1)$-fold cover of $Y=M_{ell}$.

On the other hand, over $\bar{\mathbb{F}}_p$, up to "homeomorphism" $Y_0(p)$ is two copies of Y glued at the supersingular points, so the Euler characteristic is $2\times(-1/12) - S$ (where $S$ is the "number" of supersingular elliptic curves).

However, since $Y_0(p)$ has semi-stable reduction (and is constant at infinity) over $\mathbb{Z}_p$, the special fiber is gotten from the generic fiber by contracting a bunch of "circles" to points, one for each nodal point of the special fiber; thus our second ($\mathrm{char}$ $p$) Euler characteristic is equal to our first ($\mathrm{char}$ $0$) Euler characteristic plus $S$.  Comparing and solving for $S$ gives the formula pretty quickly.