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 $M_B$ be the moduli space of the $B$-model: it is the moduli space of complex structures on $X$. Let $M_A$ 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 $M_A$ in an abstract way as space of deformation of $2d$ $TQFTs$. $M_A$ and $M_B$ 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 $F_{g,A}$ (resp. $F_{g,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. $F_{g,A}$ is a function (in fact a section of some line bundle) over $M_A$; we denote $F_{g,A}(t, \bar{t})$ where $(t, \bar{t})$ are complex coordinates over the complex variety $M_A$. The asymptotic behavior of $F_{g,A}$ around the large volume limit of $M_A$ 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 $F_{g,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: $F_1$ 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, $F_{g,A}$ is completely determined up to an holomorphic section of a line bundle. The space $M_A$ 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 $F_{g,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 $F_{1,A}$, the finite extra necessary information is empty: $F_{1,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, $F_{g,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 $M_A$ and this is obtained by the genus $0$ mirror symmetry: $M_A$ of $X$ is $M_B$ 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 $M_B$ and the genus $0$ $B$-model (it is classical complex geometry). The problem is that there is no direct definition of $M_A$ 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 $F_{g,A}$ around the large volume limit. We can define $M_A$ indirectly by mirror symmetry but then it is not clear how to extend $F_{g,A}$ to all $M_A$. 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 $F_g$'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).
Q&A (4871)
