Applications in fields like physics, biology, and chemistry often involve partial differential equation (PDE) systems with multiple interacting components, such as gas mixtures, competing populations, or chemical reactants. These systems are modeled with nonlinear reaction-diffusion equations that may include cross-diffusion terms. Developing numerical methods for these systems is challenging due to their nonlinearity, coupled nature, the need to ensure positivity and boundedness in solutions, and the often non-symmetric or non-positive definite diffusion matrix. Motivated by the inherent entropy structure of these PDE systems, we present numerical methods based on nonlinear transformations involving entropy variables to preserve boundedness in the approximated solutions. We focus on a Local Discontinuous Galerkin (LDG) method, where auxiliary variables help reformulate the problem so that nonlinearities do not appear under differential operators and interface terms. This approach allows for a parallel evaluation of nonlinear operators, supports high spatial approximation degrees, preserves boundedness of the physical variables without the need of additional postprocessing or slope limiters, and provides a discrete version of entropy stability.