In this investigation, a continuum mechanics based bilinear shear deformable shell element is developed using the absolute nodal coordinate formulation (ANCF) for the large deformation analysis of multibody shell structures. The element consists of four nodes, each of which has the global position coordinates and the transverse gradient coordinates along the thickness introduced for describing the orientation and deformation of the cross section of the shell element. The global position field on the middle surface and the position vector gradient at a material point in the element are interpolated by bilinear polynomials. The continuum mechanics approach is used to formulate the generalized elastic forces, allowing for the consideration of nonlinear constitutive models in a straightforward manner. The element lockings exhibited in the element are eliminated using the assumed natural strain (ANS) and enhanced assumed strain (EAS) approaches. In particular, the combined ANS and EAS approach is introduced to alleviate the thickness locking arising from the erroneous transverse normal strain distribution. Several numerical examples are presented in order to demonstrate the accuracy and the rate of convergence of numerical solutions obtained by the continuum mechanics based bilinear shear deformable ANCF shell element proposed in this investigation.