In this study, Genetic Programming (GP) is explored as a data-driven approach to reconstruct eddy-resolved simulations of the North Atlantic double-gyre. The North Atlantic double-gyre is simulated by a stratified quasi-geostrophic model which is solved in an eddy-resolving turbulent regime. The statistically stationary simulations of the double-gyre model are considered to analyse interdecadal variability. The simulation results are compressed using the classic Proper Orthogonal Decomposition (POD) and Spectral Proper Orthogonal Decomposition (SPOD) to characterise spatiotemporal coherent structures in mesoscale oceanic turbulence. The time variant coefficients of the two reduced-order models are next fed into a Genetic Programming code for approximation. To this end, the parameter space of objective functions in Genetic Programming is explored to capture the key statistical properties of the original time series such as variance and auto-correlation function. Finally, the flow field is reconstructed using POD and SPOD expansions with coefficients approximated by GP and the results are compared with original simulations beyond the training window.