An MHD-RE hybrid model is developed and implemented into the MARS-F code (Liu et al 2000 Phys. Plasmas7 3681), where the runaway electrons (REs) are treated as a fluid species. The model is utilized to study the tearing mode (TM) stability in a tokamak plasma, with the numerical results compared to that of the analytic theory. The RE effect on the n = 1 TM is systematically investigated in the absence or presence of the favorable average curvature stabilization (the GGJ-effect). Without the GGJ-effect, for a force-free toroidal plasma equilibrium, REs are found to destabilize the TM and induce a finite real frequency to the mode (in a plasma with vanishing equilibrium flow). The RE contribution is found to significantly enhance the toroidal mode coupling. In the presence of the GGJ-effect, for a plasma with finite equilibrium pressure, the MARS-F hybrid model finds that REs destroy the symmetry in the TM spectrum and introduce extra branch of the TM solution due to the eigenvalue bifurcation. This happens because REs lead to coupling to a stable fluid TM mode and destabilize this mode. Despite complicated mode spectrum introduced by REs, we find that the latter play an overall destabilizing role for the TM also in the presence of the GGJ-effect. The MARS-F computed eigenvalues are well fitted with the analytic dispersion relation derived by Liu et al (2020 Phys. Plasmas27 092507) for the force-free equilibrium. Less satisfactory is the fitting result with a TM dispersion relation that includes both the GGJ-effect and the RE contribution, suggesting further improvement of the proposed dispersion relation.