We investigate holomorphic anomalies of partition functions underlying string compactifications on Calabi-Yau fourfolds with background fluxes. For elliptic fourfolds the partition functions have an alternative interpretation as elliptic genera of N=1 supersymmetric string theories in four dimensions, or as generating functions for relative Gromov-Witten invariants of fourfolds with fluxes. We derive the holomorphic anomaly equations by starting from the BCOV formalism of topological strings, and translating them into geometrical terms. The result can be recast into modular and elliptic anomaly equations. As a new feature, as compared to threefolds, we find an extra contribution which is given by a gravitational descendant invariant. This leads to linear terms in the anomaly equations, which support an algebra of derivatives mapping between partition functions of the various flux sectors. These geometric features are mirrored by certain properties of quasi-Jacobi forms. We also offer an interpretation of the physics from the viewpoint of the worldsheet theory.