In this paper, we address the design of message-passing receiver for massive multiple-input multiple-output orthogonal frequency division multiplex (MIMO-OFDM) systems. With the aid of the central limit argument and Taylor-series approximation, a computationally efficient receiver that performs joint channel estimation and decoding is devised by the framework of expectation propagation. In particular, the local belief defined at the channel transition function is expanded up to the second order with Wirtinger calculus, to transform the messages sent by the channel transition function to a tractable form. As a result, the channel impulse response between each pair of antennas is estimated by Gaussian message passing. In addition, a variational expectation-maximization-based method is derived to learn the channel power-delay profiles. The proposed scheme is assessed in 3D massive MIMO-OFDM systems with spatially correlated channels, and the empirical results corroborate its superiority in terms of performance and complexity.