We study the reaction e+e-→π+π-π0ηc for the first time using data samples collected with the BESIII detector at center-of-mass energies s=4.226, 4.258, 4.358, 4.416, and 4.600 GeV. Evidence of this process is found and the Born cross section σB(e+e-→π+π-π0ηc), excluding e+e-→ωηc and ηηc, is measured to be (46-11+12±10) pb at s=4.226 GeV. Evidence for the decay Zc(3900)±→ρ±ηc is reported at s=4.226 GeV with a significance of 3.9σ, including systematic uncertainties, and the Born cross section times branching fraction σB(e+e-→π Zc(3900)±)×B(Zc(3900)±→ρ±ηc) is measured to be (48±11±11) pb, which indicates that e+e-→π Zc(3900)±→πρ±ηc dominates the e+e-→π+π-π0ηc process. The Zc(3900)±→ρ±ηc signal is not significant at the other center-of-mass energies and the corresponding upper limits are determined. In addition, no significant signal is observed in a search for Zc(4020)±→ρ±ηc with the same data samples. The ratios RZc(3900)=B(Zc(3900)±→ρ±ηc)/B(Zc(3900)±→π±J/ψ) and RZc(4020)=B(Zc(4020)±→ρ±ηc)/B(Zc(4020)±→π±hc) are obtained and compared with different theoretical interpretations of the Zc(3900)± and Zc(4020)±. © 2019 authors. Published by the American Physical Society.