This paper presents the development of accurate turbulence closures for low-pressure turbine (LPT) wake mixing prediction by integrating a machine-learning approach based on gene expression programming (GEP), with Reynolds-averaged Navier-Stokes (RANS)-based computational fluid dynamics (CFD). In order to further improve the performance and robustness of GEP-based data-driven closures, the fitness of models is evaluated by running RANS calculations in an integrated way, instead of an algebraic function. Using a canonical turbine wake with inlet conditions prescribed based on high-fidelity data of the T106A cascade, we demonstrate that the "CFD-driven"machine-learning approach produces physically correct non-linear turbulence closures, i.e., predict the right down-stream wake development and maintain an accurate peak wake loss throughout the domain. We then extend our analysis to full turbine blade cases and show that the model development is sensitive to the training region due to the presence of deterministic unsteadiness in the near-wake. Models developed including the near-wake have artificially large diffusion coefficients to over-compensate for the vortex shedding steady RANS cannot capture. In contrast, excluding the near-wake in the model development produces the correct physical model behavior, but predictive accuracy in the near-wake remains unsatisfactory. This can be remedied by using the physically consistent models in unsteady RANS. Overall, the "CFD-driven"models were found to be robust and capture the correct physical wake mixing behavior across different LPT operating conditions and airfoils such as T106C and PakB. © 2021 by ASME.