We consider the short-distance expansion of the product of two gluon field strength tensors connected by a straight-line-ordered Wilson line. The vacuum expectation value of this nonlocal operator is a common object in studies of the QCD vacuum structure, whereas its nucleon expectation value is known as the gluon quasi-parton distribution and is receiving a lot of attention as a tool to extract gluon distribution functions from lattice calculations. Extending our previous study [1], we calculate the three-loop coefficient functions of the scalar operators in the operator product expansion up to dimension four. As a by-product, the three-loop anomalous dimension of the nonlocal two-gluon operator is obtained as well.