Multidisciplinary design optimization (MDO) has become a popular tool when dealing with complex systems, since it allows to account for the interaction between the participating disciplines directly in the optimization framework. However, the use of high-fidelity solvers, based on discretized partial differential equations (PDEs), in the MDO framework remains a difficult task due to the inherent computational cost. We address this problem by replacing each disciplinary solver involved in the system by a dedicated non-intrusive surrogate approximation, obtained using Gaussian Processes. Due to the dimension of the coupling variable space (often represented by displacement and pressure fields), a dimension reduction step is proposed before the construction of the disciplinary surrogates. Global Proper Orthogonal Decomposition (POD) bases are used to reduce the dimension of the subsystems whose behavior can be described using only a few modes. For more complex behaviors, the interpolation of a set of pointwise local POD bases is preferred as dimension reduction strategy. The observed reduction in computational cost for an engineering test case using high-fidelity disciplinary solvers confirms the interest of the proposed approach. The number of disciplinary solver calls is reduced thanks to the combined dimension reduction and interpolation strategy. Comparison with the results obtained when the true disciplinary solvers are used shows that the relative error remains within acceptable bounds.