This paper presents a novel model order reduction technique for 3D flexible multibody systems featuring nonlinear elastic behavior. We adopt the mean-axis floating frame approach in combination with an enhanced Rubin substructuring technique for the construction of the reduction basis. The standard Rubin reduction basis is augmented with the modal derivatives of both free-interface vibration modes and attachment modes to consider the bending–stretching coupling effects for each flexible body. The mean-axis frame generally yields relative displacements and rotations of smaller magnitude when compared to the one obtained by the nodal-fixed floating frame. This positively impacts the accuracy of the reduction basis. Also, when equipped with modal derivatives, the Rubin method better considers the geometric nonlinearities than the Craig–Bampton method, as it comprises vibration modes and modal derivatives featuring free motion of the interface. The nonlinear coupling between free-interface modes and attachment modes is also considered. Numerical tests confirm that the proposed method is more accurate than Craig–Bampton’s, a nodal fixed floating frame counterpart originally proposed in Wu and Tiso (Multibody Syst. Dyn. 36(4): 405–425, [2016]), and produces significant speed-ups. However, the offline cost is increased because the mean-axis formulation produces operators with decreased sparsity patterns.