We propose a non-intrusive reduced basis (RB) method for parametrized nonlinear partial differential equations (PDEs) that leverages models of different accuracy. From a collection of low-fidelity (LF) snapshots, parameter locations are extracted for the evaluations of high-fidelity (HF) snapshots to recover a reduced basis. Multi-fidelity Gaussian process regression (GPR) is employed to approximate the combination coefficients of the reduced basis. LF data is assimilated either via projection onto an LF basis or via an interpolation approach inspired by bifidelity reconstruction. The correlation between HF and LF data is modeled with hyperparameters whose values are automatically determined in the regression step. The proposed methods not only leverage the assimilated LF data to reduce the cost of the offline phase, but also allow for a fast evaluation during the online stage, independent of the computational cost of neither the low-nor the high-fidelity solution. Numerical studies demonstrate the effectiveness of the proposed approach on manufactured examples and problems in nonlinear structural mechanics. Clear benefits of using lower resolution models rather than reduced physics models are observed in both the basis selection and the regression step. An active learning scheme is used for additional snapshot selection at locations with high error. The speed-up in the online evaluation and the high accuracy of extracted quantities of interest makes the multifidelity RB method a powerful tool for outer-loop applications in engineering, as exemplified in uncertainty quantification. (C) 2020 Elsevier B.V. All rights reserved.