In a recent series of papers we have shown how the eikonal/geometrical optics approximation can be used to calculate analytically the fundamental quasinormal mode frequencies associated with coupled systems of wave equations, which arise, for instance, in the study of perturbations of black holes in gravity theories beyond General Relativity. As a continuation to this series, we here focus on the quasinormal modes of nonrotating black holes in scalar Gauss-Bonnet gravity assuming a small-coupling expansion. We show that the axial perturbations are purely tensorial and are described by a modified Regge-Wheeler equation, while the polar perturbations are of mixed scalar-tensor character and are described by a system of two coupled wave equations. When applied to these equations, the eikonal machinery leads to axial modes that deviate from the general relativistic results at quadratic order in the Gauss-Bonnet coupling constant. We show that this result is in agreement with an analysis of unstable circular null orbits around blackholes in this theory, allowing us to establish the geometrical optics-null geodesic correspondence for the axial modes. For the polar modes the small-coupling approximation forces us to consider the ordering between eikonal and small-coupling perturbative parameters; one of which we show, by explicit comparison against numerical data, yields the correct identification of the quasinormal modes of the scalar-tensor coupled system of wave equations. These corrections lift the general relativistic degeneracy between scalar and tensorial eikonal quasinormal modes at quadratic order in Gauss-Bonnet coupling in a way reminiscent of the Zeeman effect. In general, our analytic, eikonal quasinormal mode frequencies (normalized to the General Relativity ones) agree with numerical results with an error of $sim 10%$ in the regime of small coupling constant. (abridged)