Developing reduced-order models for turbulent flows, which contain dynamics over a wide range of scales, is an extremely challenging problem. In statistical mechanics, the Mori-Zwanzig (MZ) formalism provides a mathematically formal procedure for constructing reduced-order representations of high-dimensional dynamical systems, where the effect due to the unresolved dynamics are captured in the memory kernel and orthogonal dynamics. Turbulence models based on MZ formalism have been scarce due to the limited knowledge of the MZ operators, which originates from the difficulty in deriving MZ kernels for complex nonlinear dynamical systems. In this work, we apply a recently developed data-driven learning algorithm, which is based on Koopmans description of dynamical systems and Moris linear projection operator, on a set of fully-resolved isotropic turbulence datasets to extract the Mori-Zwanzig operators. With data augmentation using known turbulence symmetries, the extracted Markov term, memory kernel, and orthogonal dynamics are statistically converged and the Generalized Fluctuation-Dissipation Relation can be verified. The properties of the memory kernel and orthogonal dynamics, and their dependence on the choices of observables are investigated to address the modeling assumptions that are commonly used in MZ-based models. A series of numerical experiments are then constructed using the extracted kernels to evaluate the memory effects on predictions. Results show that the prediction errors are strongly affected by the choice of observables and can be further reduced by including the past history of the observables in the memory kernel.