In general, in the approach to shadow tomography introduced in that paper, the "classical shadows" are effectively operators which have the same dimensionality as the states to be estimated (you can see them as estimators for the state itself). This means that if you think about the many-qubit scenario, you won't even be able to save this in your computer as they'll be $2^n\times 2^n$ dimensional matrices.
The stabilizer formalism allows you to represent efficiently these shadows. Rather than writing out the entire operators, you can represent them efficiently via their stabilizers. I think the main reason you can do this is that in at least the main cases studied in the paper, if I remember correctly, the shadows are proportional to the state associated to the observed outcome, so representing efficiently the state is the same as representing efficiently the shadows.
See e.g. What is a stabilizer state? for a concise explanations of how stabilizers work in general.
Let me clarify a little bit the above paragraph.
- At the measurement stage, the gist of the protocol is that you perform projective measurements in random directions. You can see this as equivalent to performing a POVM with infinitely many outcomes, with each outcome corresponding to a pure state. When this protocol is implemented via Clifford circuits, the pure states corresponding to the different outcomes can all be described in the stabilizer formalism.
- At the estimation stage, you take the different outcomes obtained from the measurement, and convert each outcome into the corresponding "classical shadow", which is an operator that at least in the cases considered in the paper is mostly proportional to the state corresponding to that outcome. These shadows are by construction such that they reproduce the correct state in expectation value, and you can show that they estimate expectation values of observables and other stuff pretty well.