Proving that $\sigma_7(n) = \sigma_3(n) + 120 \sum_{m=1}^{n-1} \sigma_3(m)\sigma_3(n-m)$ without using modular forms?

Question

This problem appears as a (starred!) exercise in D. Zagier's notes on modular forms. I have to admit that I have no idea how to do it.

Here, $\sigma_k(n) =\sum_{d\mid n} d^k$, as usual.

This identity is traditionally obtained by using the fact that the space $M_8(\text{SL}_2(\mathbf Z))$ of modular forms of weight $8$ and level $1$ is $1$-dimensional, and contains both $E_4^2$ and $E_8$ ($E_k$ the Eisenstein series of weight $k$). Using this, it's a piece of cake (simply a matter of comparing coefficients).

Without using modular forms, though, I am stumped.