Graphene bilayers display peculiar electronic and mechanical characteristics associated with their two-dimensional character and relative disposition of the sheets. Here, we study nuclear quantum effects in graphene bilayers by using path-integral molecular dynamics simulations, which allow us to consider quantization of vibrational modes and study the effect of anharmonicity on physical variables. Finite-temperature properties are analyzed in the range from 12 to 2000 K. Our results for graphene bilayers are compared with those found for graphene monolayers and graphite. Nuclear quantum effects turn out to be appreciable in the layer area and interlayer distance at finite temperatures. Differences in the behavior of in-plane and real areas of the graphene sheets are discussed. The interlayer spacing has a zero-point expansion of 1.5 × 10-2 Å with respect to the classical minimum. The compressibility of graphene bilayers in the out-of-plane direction is found to be similar to that of graphite at low temperatures and increases faster as the temperature is raised. The low-temperature compressibility increases by 6% due to zero-point motion. Special emphasis is placed on atomic vibrations in the out-of-plane direction. Quantum effects are present in these vibrational modes, but classical thermal motion becomes dominant over quantum delocalization for large system size. The significance of anharmonicities in this atomic motion is estimated by comparing with a harmonic approximation for the vibrational modes in graphene bilayers.