In this paper, a nonlinear time-varying dynamic model of a drivetrain composed of a spiral bevel gear pair, shafts and bearings is developed. Gear shafts are modeled by utilizing Timoshenko beam finite elements, and the mesh model of a spiral bevel gear pair is used to couple them. The dynamic model includes the flexibilities of shaft bearings as well. Gear backlash and time variation of mesh stiffness are incorporated into the dynamic model. Clearance nonlinearity of bearings is assumed to be negligible, which is valid for preloaded rolling element bearings. Furthermore, stiffness fluctuations of bearings are disregarded. Multi-term harmonic balance method (HBM) is applied on the system of nonlinear differential equations in order to obtain a system of nonlinear algebraic equations. Utilizing receptance method, system of nonlinear algebraic equations is grouped in nonlinear and linear sets of algebraic equations where the nonlinear set can be solved alone decreasing the number of equations to be solved significantly. This reduces the computational effort drastically which makes it possible to use finite element models for gear shafts. In the calculation of Fourier coefficients, continuous-time Fourier transform as opposed to the gear dynamics studies that utilize discrete Fourier Transform is used. Thus, convergence problems that arise when the number of nonlinear DOFs is large are avoided. Moreover, analytical integration is employed for the calculation of Fourier coefficients rather than numerical integration in order to further reduce the computational time required. Nonlinear algebraic equations obtained are solved by utilizing Newton's method with arc-length continuation. Direct numerical integration is employed to verify the solutions obtained by HBM. Several case studies are carried out, and the influence of backlash amount, fluctuation of gear mesh stiffness and variation of bearing stiffness are investigated. In addition to these, the response of the coupled gear system model is compared with that of gear torsional model in order to study the influence of the coupling on dynamics of the system.