The increasing demands for travelling comfort and reduction of carbon dioxide emissions have been considered substantially in the stage of conceptual aircraft design. However, the design of a modern aircraft is a multidisciplinary process, which requires the coordination of information from several specific disciplines, such as structures, aerodynamics, control, etc. To address this problem with adequate accuracy, the multidisciplinary analysis and optimization (MAO) method is usually applied as a systematic and robust approach to solve such complex design issues arising from industries. Since MAO method is tedious and computationally expensive, genetic programming (GP)-based metamodeling techniques incorporating MAO are proposed as an effective approach to minimize the wing stiffness of a large aircraft subject to aerodynamic, aeroelastic and stability constraints in the conceptual design phase. Based on the linear small-disturbance theory, the state-space equation is employed for stability analysis. In the process of multidisciplinary analysis, aeroelastic response simulations are performed using Nastran. To construct metamodels representing the responses of the interests with high accuracy as well as less computational burden, optimal Latin hypercube design of experiments (DoE) is applied to determine the optimized distribution of sampling points. Following that, parametric optimization is carried out on metamodels to obtain the optimal wing geometry shape, elastic axis positions and stiffness distribution, and then the solution is verified by finite element simulations. Finally, the superiority of the GP-based metamodel technique over genetic algorithm is demonstrated by multidisciplinary design optimization of a representative beam-frame wing structure in terms of accuracy and efficiency. The results also show that GP metamodel-based strategy for solving MAO problems can provide valuable insights to tailoring parameters for the effective design of a large aircraft in the conceptual phase.