Multiscale topology optimization (MSTO) is a numerical design approach to optimally distribute material within coupled design domains at multiple length scales. Due to the substantial computational cost of performing topology optimization at multiple scales, MSTO methods often feature subroutines such as homogenization of parameterized unit cells and inverse homogenization of periodic microstructures. Parameterized unit cells are of great practical use, but limit the design to a pre-selected cell shape. On the other hand, inverse homogenization provide a physical representation of an optimal periodic microstructure at every discrete location, but do not necessarily embody a manufacturable structure. To address these limitations, this paper introduces a Gaussian process regression model-assisted MSTO method that features the optimal distribution of material at the macroscale and topology optimization of a manufacturable microscale structure. In the proposed approach, a macroscale optimization problem is solved using a gradient-based optimizer The design variables are defined as the homogenized stiffness tensors of the microscale topologies. As such, analytical sensitivity is not possible so the sensitivity coefficients are approximated using finite differences after each microscale topology is optimized. The computational cost of optimizing each microstructure is dramatically reduced by using Gaussian process regression models to approximate the homogenized stiffness tensor. The capability of the proposed MSTO method is demonstrated with two three-dimensional numerical examples. The correlation of the Gaussian process regression models are presented along with the final multiscale topologies for the two examples: a cantilever beam and a 3-point bending beam.