Reynolds-averaged Navier–Stokes (RANS) models remain the backbone of industrial computational fluid dynamics due to their robustness and low computational cost, yet their turbulence closures rely on structural assumptions that limit accuracy in separated, anisotropic, and non-equilibrium flows. Data-driven approaches offer a route to reduce these model-form errors, but many existing methods rely on black-box models or training targets that are inconsistent with the RANS operator, leading to limited interpretability and fragile a posteriori performance. This work presents a data-driven turbulence-model correction framework based on symbolic regression. Instead of learning Reynolds stresses directly from high-fidelity data, model discrepancies are represented as an explicit corrective forcing term in the RANS momentum equations. The forcing is extracted through a variational formulation that ensures consistency with the discretized RANS equations, turbulence variables, and mesh resolution. Local, dimensionless flow features constructed from normalized strain, rotation, pressure-gradient, and turbulence scales are used as inputs. A deep neural network is first employed as a reference model to establish achievable accuracy and to identify redundant features. Building on this reduced feature set, symbolic regression is used to discover compact analytical expressions for the corrective forcing. The approach is demonstrated on canonical turbulent flow benchmarks, including the periodic hill configuration. The resulting models are stable and expressed in closed analytical form. Ongoing and future work focuses on the direct integration of the discovered corrections into OpenFOAM-based RANS solvers, enabling fully embedded simulations without external inference, and on assessing robustness across Reynolds numbers and flow configurations.