In this paper, we study time-varying channel estimation and optimal pilot design for multiple-input multiple-output two-way multi-relay systems in the sense of minimizing the total estimation mean square error (MSE) under the power constraints at two source nodes and multiple relays. A particularly challenging issue in the optimal pilot design is to completely eliminate the inter-link interference while most efficiently using the channel resource. To address this challenging issue, we consider a pilot design approach based on the phase rotation technique at multiple relays. First, we establish various optimality conditions for the pilot design in the sense of the minimum total MSE. Then, we propose the optimal pilot scheme by designing the pilot signals, phase rotations, and pilot symbol positions to satisfy the established optimality conditions. Finally, we propose the asymptotically optimal pilot scheme in the high signal-to-noise ratio (SNR) regime with low complexity. Simulation results show that the proposed schemes significantly outperform the existing schemes in terms of the channel estimation and the bit error rate, and the proposed asymptotically optimal scheme provides the near-optimal performance even in the low to moderate SNR range.