Theorem. If $\alpha,\beta>0$, $\alpha+\beta<n$, then $I_\alpha(I_\beta\varphi)=I_{\alpha+\beta}\varphi$ for $\varphi\in\mathscr S_n$.
Lemma. If $\alpha,\beta>0$, $\alpha+\beta<n$, then there is a constant $C_0=C_0(\alpha,\beta,n)$ such that $$ \int_{\mathbb{R}^n}\frac{dy}{|x-y|^{n-\alpha}|y|^{n-\beta}}= \frac{C_0}{|x|^{n-(\alpha+\beta)}}\, . $$
Proof. First you show that the integral is finite for every $x\neq 0$. By rotational symmetry, the integral on the left hand side depends on $|x|$ only. Denoting its value by $f(|x|)$ a simple change of variables (by scaling) show that $f(|x|)=|x|^{\alpha+\beta-n}f(1)$ and the result follows. $\Box$
Proof of the theorem. The lemma and the Fubini theorem easily implies that $$ I_\alpha(I_\beta\varphi)(x)= \frac{C_0\gamma(\alpha+\beta)}{\gamma(\alpha)\gamma(\beta)} I_{\alpha+\beta}\varphi(x). $$ The only problem is to show that the constant is actually equal $1$.
To prove this it suffices to verify that $I_\alpha(I_\beta\varphi)= I_{\alpha+\beta}\varphi$ for just one non-zero function $\varphi$. To this end let $\varphi\in\mathscr S_n$ be such that $\hat{\varphi}=0$ in a neighborhood of $0$. Then $$ I_\alpha(I_\beta\varphi)= I_\alpha\Big(\Big(\underbrace{(4\pi^2|\xi|^2)^{-\beta/2}\hat{\varphi}}_{\in\mathscr S_n}\Big)^\vee\Big) = \left((4\pi^2|\xi|^2)^{-\alpha/2}(4\pi^2|\xi|^2)^{-\beta/2}\hat{\varphi}\right)^\vee= I_{\alpha+\beta}\varphi. $$ The proof is complete. $\Box$