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.