Vesselin has already given an excellent answer to your first question, so let me just elaborate a little bit on the relation
$$
28\zeta(3,9)+150\zeta(5,7)+168\zeta(7,5)=\frac{5197}{691}\zeta(12)\tag{1}
$$
and its connection with modular forms for $\operatorname{SL}_2(\mathbb{Z})$. First of all, it is a rather exotic and surprising relation (I will qualify this statement below). In a nutshell, the coefficients $28,150,168$ on the left-hand side of (1) come from period integrals of the cusp form $\Delta$, while the denominator $691$ of the right-hand side comes from the constant term of the Eisenstein series $G_{12}$.
To any cuspidal Hecke eigenform $f$ of weight $k$, one can associate an even (i.e. invariant under $Y\mapsto -Y$) homogeneous polynomial $r_f(X,Y)^{+} \in K_f[X,Y]_{k-2}$ of degree $k-2$, where $K_f$ is the number field generated by the Fourier coefficients of $f$. Explicitly,
$$
r_f(X,Y)^{+}=\frac{\operatorname{Re}\left( (2\pi i)^{k-1}\int_0^{i\infty}f(z)(X-zY)^{k-2}dz \right)}{\omega_f^{+}}
$$
where the integral is along the geodesic path from $0$ to $i\infty$ on the upper half-plane, and $\omega_f^{+} \in \mathbb{R}$ is a suitable ''period''. This construction is part of the very well-known Eichler--Shimura--Manin theory, expositions of which can be found in many places; see e.g. Lang's Introduction to Modular Forms for an introduction, or Brown's Multiple Modular Values and the relative completion of the fundamental group of $M_{1,1}$, Section 7 for a more condensed presentation.
Now in weight $12$ there is a unique cuspidal Hecke eigenform, namely Ramanujan's $\Delta$, and one has (up to scaling by an element of $\mathbb{Q}^{\times}$)
$$
r_{\Delta}(X,Y)^+=\frac{36}{691}(X^{10}-Y^{10})-X^2Y^2(X^2-Y^2)^3.
$$
The essential part is the second summand $p(X,Y):=X^2Y^2(X^2-Y^2)^3$, the first one is a coboundary in the sense of group cohomology for $\operatorname{SL}_2(\mathbb{Z})$, and hence can be ignored (for the time being, but see below). Now Gangl--Kaneko--Zagier's result (Theorem 3 in their paper) says in essence that if numbers $q_{r,s}$ are defined by $p(X+Y,Y)=\sum_{r,s \geq 1}\binom{r+s-2}{r-1}q_{r,s}X^{r-1}Y^{s-1}$, then
$$
\sum_{r+s=12, \, r,s \, odd}q_{r,s}\zeta(r,s) \equiv 0 \mod \zeta(12)
$$
which, after rescaling, yields (1) modulo $\zeta(12)$. The upshot is that the coefficients $28,150,168$ come from integrals involving the cusp form $\Delta$, and therefore the left-hand side of (1) has a rather clean, arithmetic-geometric interpretation.
To the best of my knowledge, no such interpretation is known for the right-hand prefactor $\frac{5197}{691}$ (in Gangl--Kaneko--Zagier's original approach, they need to compute some kind of beta-integral). However, the occurrence of the denominator $691$ seems to be an artifact of the well-known congruence $G_{12} \equiv \Delta \mod 691$ on the level of Fourier expansions, where $G_{12}$ is the Hecke Eisenstein series of weight $12$ (basically because the numerator of the constant term of G_{12} is divisible by 691). The precise relation should be explained by recent work of Hain--Matsumoto on Universal Mixed Elliptic Motives (in particular, the previously neglected term $\frac{36}{691}(X^{10}-Y^{10})$ resurfaces as the dual of the cocycle of $G_{12}$). As far as I can see, this does not give an explanation about the numerator $5197$, though.
Finally, I would like to point out that there is another connection of (1) to cusp forms, which is not (yet?) arithmetic-geometric but in my opinion very nice: namely, via the theory of double Eisenstein series, which are to usual Eisenstein series what double zeta values are to single zeta values. Since my answer is already much too long, I refer to the original paper of Gangl--Kaneko--Zagier, Section 7 for the details. Some of this has also been generalized by Henrik Bachmann in his PhD thesis (see in particular the bottom of p.34 for details on how to obtain (1) via this theory).
So why is the relation (1) surprising? It is more generally believed that all algebraic relations between multiple zeta values are implied by the (regularized) double shuffle equations, which are a class of quadratic relations between multiple zeta values. Specialized to the depth two case, they read
$$
\zeta(r)\zeta(s)=\zeta(r,s)+\zeta(s,r)+\zeta(r+s)=\sum_{m+n=r+s}\left( \binom{m-1}{r-1}+\binom{m-1}{s-1}\right)\zeta(m,n). \tag{2}
$$
One of the main results of Gangl--Kaneko--Zagier is that (1) cannot be deduced from (2), and one needs to resort to the double shuffle equations in higher depth. If one were able to deduce from (2) that $28\zeta(4)\zeta(8)+\frac{95}{3}\zeta(6)^2 \in \mathbb{Q}\zeta(12)$ (which is of course a consequence of Euler's theorem on even zeta values), this ''depth-defect'' would go away, alas one cannot. This is one of the very perplexing features of the depth filtration on multiple zeta values, and that it can be resolved using modular forms points at some deep relation between the two objects, which we do not completely understand yet.
