Due to the lack of knowledge in terms of their flexibility and deformation, spline joints are typically assumed to be rigid in dynamic models of gearboxes, transmissions, and drivetrains. As various dynamic phenomena are associated with the stiffness of a spline joint, any high-fidelity dynamic model of drivetrains must properly capture the stiffness of spline joints. In this study, a general analytical stiffness formulation for spline joints is proposed based on a semi-analytical spline load distribution model. This formulation defines a fully populated stiffness matrix of a spline joint including radial, tilting, and torsional stiffness values as well as off-diagonal coupling terms. A blockwise inversion method is proposed and implemented with this analytical formulation to reduce computational time required. At the end, a detailed parametric study is presented to demonstrate the sensitivity of the spline stiffness matrix to torque level, tooth modifications, misalignments, and tooth indexing errors.