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 char. 0 and in char. p, then see how they must relate. Here we go:
Over C, and hence over Q_p-bar as well, the Euler characteristic of Y_0(p) is (p+1)*(-1/12), since Y_0(p) is a (p+1)-fold cover of Y=M_{ell}.
On the other hand, over F_p-bar, up to "homeomorphism" Y_0(p) is two copies of Y glued at the supersingular points, so the Euler characteristic is 2*(-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 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 (char p) Euler characteristic is equal to our first (char 0) Euler characteristic plus S. Comparing and solving for S gives the formula pretty quickly.