Using a data sample of 921.9 fb$^{-1}$ collected with the Belle detector, we study the process of $e^+e^-to D^+_sD_{s1}(2536)^-+c.c.$ via initial-state radiation. We report the first observation of a vector charmoniumlike state decaying to $D^+_sD_{s1}(2536)^-+c.c.$ with a significance of 5.9$sigma$, including the systematic uncertainties. The measured mass and width are $(4625.9^{+6.2}_{-6.0}({rm stat.})pm0.4({rm syst.}))~{rm MeV}/c^{2}$ and $(49.8^{+13.9}_{-11.5}({rm stat.})pm4.0({rm syst.}))~{rm MeV}$, respectively. The product of the $e^+e^-to D^+_sD_{s1}(2536)^-+c.c.$ cross section and the branching fraction of $D_{s1}(2536)^-to{bar D}^{*0}K^-$ is measured from the $D_s bar{D}_{s1}(2536)$ threshold to 5.59~GeV.