We apply a multiscale modeling approach to study lattice reconstruction in marginally twisted bilayers of transition metal dichalcogenides (TMD). For this, we develop DFT-parametrized interpolation formulae for interlayer adhesion energies of MoSe\(_2\), WSe\(_2\), MoS\(_2\), and WS\(_2\), combine those with elasticity theory, and analyze the bilayer lattice relaxation into mesoscale domain structures. Paying particular attention to the inversion asymmetry of TMD monolayers, we show that 3R and 2H stacking domains, separated by a network of dislocations develop for twist angles \(\theta^{\circ}<\theta^{\circ}_P\sim 2.5^{\circ}\) and \(\theta^{\circ}<\theta^{\circ}_{AP}\sim 1^{\circ}\) for, respectively, bilayers with parallel (P) and antiparallel (AP) orientation of the monolayer unit cells and suggest how the domain structures would manifest itself in local probe scanning of marginally twisted P- and AP-bilayers.