Gravitons couple to energy-momentum tensor $T_{\mu \nu}$, not mass (which captures only the $00$ component of $T_{\mu \nu}$). That's why photons (which are also massless, like gravitons) interact with gravitons leading to bending of light (Gravitational lensing).
The idea is that a free graviton field has its own $T_{\mu \nu}$, which must be added to the matter source's $T_{\mu \nu}$ to ensure that the conservation equation/Bianchi identity holds. But note that this free graviton's $T_{\mu \nu}$ is quadratic in the graviton field (say $h_{\mu \nu}$). This $O(h^2)$ term in the equations of motion can only come from an $O(h^3)$ term in the Lagrangian. This means that we have to update our free graviton theory to one that has a cubic interaction. But then this $O(h^3)$ term in the Lagrangian produces an $O(h^3)$ energy-momentum tensor, which must again be added to the previously-obtained sum of energy-momentum tensors in the equations of motion. This $O(h^3)$ term in the equations of motion can only come from an $O(h^4)$ term in the Lagrangian. This leads to a quartic interaction for the graviton field. And so on, until we sum up the entire series to obtain non-linear Einstein's theory with $n$-point graviton-graviton interactions for every $n$ from $n = 3$ until infinity. This was discussed several decades ago, and is nicely summarized in Deser's 1969 paper.