I begin with a physical summary.
Let X be a Calabi-Yau 3-fold. The supersymmetric 2d sigma model of target X defines a N=(2,2) 2d SCFT (SuperConformal Field Theory).There are essentially two ways to topologically twist this theory, the A and B-twists giving rise to the A and B-models which are 2d TQFTs (Topological Quantum Field Theory). The A-model depends only on the Kähler structure on X whereas the B-model depends only on the complex structure on X. Let MB be the moduli space of the B-model: it is the moduli space of complex structures on X. Let MA be the moduli space of the A-model. In the large volume limit of X, it coincides with the complexified Kähler cone of X. But for small Kähler classes, the sigma model description breaks down but we can extend the definition of MA in an abstract way as space of deformation of 2d TQFTs. MA and MB are special Kähler manifolds.
Coupling of the A- (resp. B-) model with 2d topological gravity defines A-(resp. B-) models of the topological string theory. Coupling to 2d gravity in genus g means essentially to integrate over the moduli space of genus g Riemann surface the correlation functions of the 2d TQFT on genus g Riemann surfaces. In particular, for each genus g, one can define Fg,A (resp. Fg,B) a genus g free energy of the A- (resp. B-) model of the topological string theory.
Let us consider for example the A-model case. Fg,A is a function (in fact a section of some line bundle) over MA; we denote Fg,A(t,ˉt) where (t,ˉt) are complex coordinates over the complex variety MA. The asymptotic behavior of Fg,A around the large volume limit of MA is related to the genus g Gromov-Witten invariants of X (i.e. approximatively the counting of genus g curves in X). The important point is that Fg,A is almost holomorphic (i.e. only function of t) but the origin of the non-holomorphy is exactly known (the non-compacity of the moduli space of genus g Riemann surfaces i.e. the boundaries of the compactification of the moduli space) and there is a formula giving a way to determine this non-holomorphic part: the holomorphic anomaly equation (first derived in genus 1 by the paper cited in the question and then for any genus in the follow-up http://arxiv.org/abs/hepth/9309140 ). (Remark: F1 is somewhat special and is in fact almost the sum of something holomorphic and of something antiholomorphic, the holomorphic anomaly equation determining the remainder).
So, thanks to the holomorphic equation, Fg,A is completely determined up to an holomorphic section of a line bundle. The space MA is in general non-compact but if we fix boundaries conditions, the set of holomorphic sections of the line bundle satisfying these conditions is finite dimensional. This means that to completely determine Fg,A, we have to know the boundary conditions plus some finite extra information. The study of boundaries condition is a local problem and so can be done. For the quintic, there are three special points to look at: the "infinity", the "orbifold point" and the "conifold point". It happens that for F1,A, the finite extra necessary information is empty: F1,A is completely determined. For higher genus, there is some necessary extra information and there is no known systematic way to proceed. The state of the art is that for the quintic, Fg,A is known for g=0,1,...,51: http://arxiv.org/abs/hep-th/0612125 , but not for g > 52.
To apply the preceding strategy, we have to know MA and this is obtained by the genus 0 mirror symmetry: MA of X is MB of the mirror of X. But the rest of the argument does not use mirror symmetry: it uses the holomorphic anomaly equation, true for both the A- and B-models, and then a global argument which is essentially basic complex analysis: a holomorphic function is more or less determined by its singularities and boundary behaviors.
The OP asks for a mathematical explanation. The problem is that few of the preceding steps are now turned into rigorous mathematics. It is probably possible to give an axiomatic definition of a 2d SCFT but it is not rigorously known (as far as I know) how to associate a 2d SCFT to a Calabi-Yau X. There is a rather well known axiomatic definition of 2d TQFT and it is possible to rigorously construct MB and the genus 0 B-model (it is classical complex geometry). The problem is that there is no direct definition of MA away from the large volume limit. The fact that we know how to define Gromov-Witten invariants means that we know how to define the genus g A-model Fg,A around the large volume limit. We can define MA indirectly by mirror symmetry but then it is not clear how to extend Fg,A to all MA. More generally, in the math litterature, mirror symmetry is often considered as a limit statement and the global geometry of the moduli spaces, which is the key of the physical argument, is not very studied (as far as I know). Nevertheless, the almost holomorphy of the Fg's has given birth to a new mathematical notion: the quasimodular forms (see Zagier's work for example).
The physical prediction for genus 1 Gromov-Witten invariants of the quintic has been mathematically proved in http://arxiv.org/abs/0705.2397 . This is not a "conceptual" proof but a direct computation (use the fact that the quintic is naturally embedded in a more "symmetric" space: the projective space, endowed with a natural torus action, and then use toric localization).