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?

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.

Actually there is a purely elementary proof of this identity due to a lemma of Liouville in his work (1858) (Liouville, J. Sur quelques formules générales qui peuvent être utiles dans la théorie des nombres). It would be helpful to understand his series of articles there (I found the third one very good) or to immitate the reproduction of proof in http://www.numdam.org/article/JTNB_2006__18_1_223_0.pdf to prove the following elementary fact on arithmetic functions (or "order theory"): If $$f:\ \mathbb{Z}^4\rightarrow\mathbb{C}$$ is a function such that $$f(a,b,x,y)-f(x,y,a,b)=f(-a,-b,x,y)-f(x,y,-a,-b)$$ for any integers $a,b,c,d$, then we have the following identity for any positive integer $n$ $$ \begin{aligned} &\sum_{ax+by=n,\\ a,b,x,y\in\mathbb{Z}_{>0}}(f(a,b,x,-y)-f(a,-b,x,y)+f(a,a-b,x+y,y)-f(a,a+b,y-x,y)+f(b-a,b,x,x+y)-f(a+b,b,x,x-y))\\ =&\sum_{d|n}\ \sum_{x\in\mathbb{Z},0<x<d}(f(o,n/d,x,d)+f(n/d,0,d,x)+f(n/d,n/d,d-x,-x)-f(x,x-d,n/d,n/d)-f(x,d,0,n/d)-f(d,x,n/d,0)). \end{aligned}$$ Now consider the function $$f(a,b,x,y)=xy^5+x^5y-2x^3 y^3$$ which is independent of $a,b$ and clearly satisfies the above condition. So we can apply the lemma above to write down the following identity. The LHS of it is $$\begin{aligned} &\sum_{ax+by=n,\\ a,b,x,y\in\mathbb{Z}_{>0}}(-2xy^5-2x^5 y+4x^3 y^3 +(x+y)y^5 + (x+y)^5y - 2(x+y)^3 y^3-(y-x)y^5-(y-x)^5y+ 2(y-x)^3 y^7+ x(x+y)^5+ x^5(x+y)-2x^3 (x+y)^3- x(x-y)^5- x^5(x-y)+ 2x^3(x-y)^3)\\ =&\sum_{ax+by=n,\\ a,b,x,y\in\mathbb{Z}_{>0}}36x^3 y^3\\ =&36\sum_{m=1}^{n-1}\sigma_3 (m)\sigma_3 (n-m). \end{aligned}$$ and the RHS is simply $$\sum_{d|n}\ \sum_{x\in\mathbb{Z},0<x<d}(xd^5+5x^2d^4-12x^3d^4+4x^4d^2+2x^5d).$$ Note that if we sum $x$ over $1$ to $d$ the result is the same, and we have the formula of sum of powers $$S_p(d)=\sum_{k=1}^{d}k^p=\frac{1}{p+1}\sum_{j=0}^{p}(-1)^j {p+1\choose j}B_j d^{p+1-j},$$ In particular one has $$\begin{aligned} S_1(d)&=\frac{d^2+d}{2}\\ S_2(d)&=\frac{2d^3+3d^2+d}{6}\\ S_3(d)&=\frac{d^4+2d^3+d^2}{4}\\ S_4(d)&=\frac{6d^5+15d^4+10d^3-d}{30}\\ S_5(d)&=\frac{2d^6+6d^5+5d^4-d^2}{12}. \end{aligned}$$ The RHS can be then simplified in terms of $S_i(d) $ ($1\leq i\leq 5$)into $$\sum_{d|n}\frac{3}{10}(d^7-d^3)=\frac{3}{10}(\sigma_7 (n)-\sigma_3 (n))$$ and thus from $LHS=RHS$ one concludes that $$\sum_{m=1}^{n-1}\sigma_3(m) \sigma_3(n-m)=\frac{1}{120}(\sigma_7(n)-\sigma_3(n)).$$ It is also interesting to note that similarly one can prove many identities arising from modular forms, more precisely from the linear dependence of elements in the algebra of modular forms of a certain weight, by using only combinatorial "A=B" tricks and "no analysis". For example, for weight $14$ one also has identities like $$ \begin{aligned} \sum_{m=1}^{n-1}\sigma_5(m) \sigma_7(n-m)&=\frac{1}{10080}(\sigma_{13}(n)+20\sigma_7(n)-21\sigma_5(n)\\ \sum_{m=1}^{n-1}\sigma_3(m) \sigma_9(n-m)&=\frac{1}{2640}(\sigma_{13}(n)-11\sigma_9(n)+10\sigma_3(m)\\ \sum_{m=1}^{n-1}\sigma(m) \sigma_{11}(n-m)&=\frac{1}{65520}(691\sigma_{13}(n)+2730(1-n)\sigma_{11}(n)-691\sigma(n)). \end{aligned}$$