In this study, the two-echelon vehicle routing problem with simultaneous pickup and delivery (2E-VRPSPD) is considered. A two-index node-based mixed integer programming (MIP) formulation is developed for the problem and then valid inequalities are used to strengthen the formulation. Moreover, several variants of the 2E-VRPSPD are introduced and the MIP formulation is adapted for these variants. To solve the problem and its variants, a matheuristic algorithm based on variable neighborhood descent algorithm with local search and mathematical programming is proposed. The performance of the proposed matheuristic is analyzed on 2E-VRPSPD and each variant of the problem using test problems derived from the literature. The experimental studies indicate that 390 out of 564 test instances up to 10 depots and 100 customers are solved to optimality for the base problem 2E-VRPSPD. Similar satisfactory results are also obtained for the other variants of the problem using same data sets.