I do not know a reference. The proof is straightforward, at least when the characteristic $p$ equals $0$ or is $>d$.
<B>Counterexample in the "wild" case.</B> When the characteristic $p$ divides one of the integers $\lambda_i$, the result is false, because of moduli of ramification. For instance, consider the $1$-parameter family of morphisms of degree $d=p+1$ from the projective line to itself, $$u_t:\mathbb{P}^1_k \to \mathbb{P}^1_k, \ \ u_t([x,y]) = [x^{p+1},xy^p+ty^{p+1}], \ t\in k^\times.$$ This is a varying family in $t$, but the ramification profile and image divisor in the target $\mathbb{P}^1_k$ is fixed. Thus, assume that $p$ equals $0$ or that $p>d$. Under this hypothesis, every finite cover $f:C\to \mathbb{P}^1_k$ is <i>tamely ramified</i>. The $k$-stack of stable maps to $\mathbb{P}^1_k$ is a proper, finitely presented Artin $k$-stack with finite diagonal. Inside this Artin $k$-stack, there is a maximal open substack that is a Deligne-Mumford $k$-stack, $$\overline{\mathcal{M}}_{g,0}(\mathbb{P}^1_k/k,d)^{\text{DM}} \subset \overline{\mathcal{M}}_{g,0}(\mathbb{P}^1_k/k,d).$$ The hypothesis implies that the open substack $\overline{\mathcal{M}}_{g,0}(\mathbb{P}^1_k/k,d)^{\text{DM}}$ equals the entire stack $\overline{\mathcal{M}}_{g,0}(\mathbb{P}^1_k/k,d)$. In particular, this open substack contains $\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)$, the locus of stable maps with smooth domain.
<B>Splitting ramification.</B> For every finite and flat morphism with Noetherian source and target, $$f:S\to T,$$ consider the coherent $\mathcal{O}_S$-module of relative differentials $\Omega_f$. This is $f$-finite, but it typically is not $f$-flat.
<B>Definition 1.</B> The $f$-<b>splitting stratification</b> is the flattening stratification of $\Omega_f$. This is a locally closed immersion that is bijective on field-valued points, $$\tau_f:T_f\hookrightarrow T,$$ such that for every morphism $u:U\to T$, there is a factorization $u=\tau_f\circ \upsilon$ for a morphism $\upsilon:U\to T_f$ if and only if the pullback of $\Omega_f$ is flat over $U$.
<B>Hypothesis 2.</B> The morphism $f$ is <b>curvilinear</b> and <b>tame</b>, i.e., every connected component of every fiber of $f$ over every geometric point of $T$, say $\text{Spec}\ k \to T$, is of the form $\text{Spec}\ k[s]/\langle s^{m+1}\rangle$ for an integer $m\geq 0$ that is prime to $p$.
<B>Notation 3.</B> The fiber product of $f$ and $\tau_f$ is denote $S_f$, and the projection morphism is denoted $\widetilde{f}:S_f\to T_f$. The pullback of $\Omega_f$ to $S_f$ is denoted $\widetilde{\Omega}_f$. The annihilator of $\widetilde{\Omega}_f$ as an ideal sheaf in $S_f$ is denote $\widetilde{\mathcal{I}}_f$. The annihilator of the ideal sheaf $\widetilde{\mathcal{I}}_f$ as an ideal sheaf in $S_f$ is denoted $\mathcal{I}_f$. The associated closed immersion is denoted $\sigma_f:S^0_f\hookrightarrow S_f$.
Under the hypotheses, $S^0_f$ is finite and étale over $T_f$, and $\sigma_f$ is bijective on field-valued points.
<B>Definition 4.</B> When $f$ is curvilinear, the commutative diagram, $$\begin{array}{ccc}S_f^0 & \hookrightarrow & S \\ \downarrow & & \downarrow \\ T_f & \hookrightarrow & T \end{array},$$ is a <b>splitting of the ramification of</b> $f$.
<B>Stable maps.</B> Denote the universal stable map as follows, $$(\pi,u):\mathcal{C}\to \mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d) \times_{\text{Spec}\ k}\mathbb{P}^1_k.$$ Denote by $\Omega$ the associated relative differentials of $(\pi,u)$ considered as a coherent sheaf on $\mathcal{C}$. Again by the hypothesis that $p>d$, the sheaf $\Omega$ is finite relative to $\pi$, i.e., the support of $\Omega$ maps finitely to $\mathcal{M}_{g,0}$ under $\pi$. In fact, $\Omega$ is $\pi$-flat and locally free of rank $b=2g+2d-2$. As a coherent sheaf on $\mathcal{C}$, the sheaf $\Omega$ is cyclic. Denote by $\mathcal{J}$ the annihilator ideal of $\Omega$ as an ideal sheaf on $\mathcal{C}$, and denote by $\mathcal{R}\subset \mathcal{C}$ the closed immersion associated to $\mathcal{J}$. Thus, $\Omega$ is the pushforward from $\mathcal{R}$ of an invertible sheaf. Denote by $\rho$ the restriction of $\pi$ to $\mathcal{R}$, $$\rho:\mathcal{R}\to \mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d).$$ This is a finite and flat morphism of degree $b=2g+2d-2$.
As in Definition 4, denote a splitting of the ramification of $\rho$ as follows,
$$
\begin{array}{ccc} \mathcal{R}^0_\rho & \xrightarrow{i} & \mathcal{R} \\
\widetilde{\rho} \downarrow & & \downarrow \rho \\
\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_\rho & \xrightarrow{\tau_\rho} & \mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d) \end{array}.
$$
This is not quite enough. The composition of $i$ with the universal stable map defines a finite morphism, $$(\widetilde{\rho},\widetilde{i}):\mathcal{R}^0_\rho\to \mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_\rho \times_{\text{Spec}\ k}\mathbb{P}^1_k.$$ Denote the image closed immersion by $$(\varpi,\iota):\mathcal{B}_\rho\hookrightarrow \mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_\rho \times_{\text{Spec}\ k}\mathbb{P}^1_k.$$ The combination of the flattening stratification of the finite morphism $\varpi$ and the splitting of the ramification of $\varpi$ gives a commutative diagram,
$$
\begin{array}
\mathcal{B}^0_{\rho,\varpi} & \xrightarrow{j} & \mathcal{B}_\rho \\
\widetilde{\varpi} \downarrow & & \downarrow \varpi \\
\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{\rho,\varpi} & \xrightarrow{\tau_{\varpi}} & \mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{\rho}
\end{array}.
$$
The composite locally closed immersion, $$\tau_{\rho,\varpi}:\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{\rho,\varpi} \to \mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d),$$ is a bijection on field-valued points. The pullback of the universal stable map via this base change has both split ramification locus in the domain curve and split branch locus in the domain curve. Thus, on each connected component of $\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{\rho,\varpi}$, there is an associated integer $n$ and an unordered $n$-tuple of partitions of $d$, say $\underline{\lambda} = \{ \lambda_1,\dots, \lambda_n\}$ with lengths $\ell_i =\ell(\lambda_i)$ such that $b$ equals the sum of the weights $$w_i = w(\lambda_i) = d -\ell(\lambda_i) = d-\ell_i.$$ Denote by $\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{n,\underline{\lambda}}$ the union of all connected components of $\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{\rho,\varpi}$ having this datum. This is the maximal open and closed such that every stable map parameterized by $\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{n,\underline{\lambda}}$ has $n$ branch points and has ramification profile $\underline{\lambda}$.
The branch locus $\mathcal{B}^0_{\rho,\varpi}$ is $\widetilde{\varpi}$-flat of relative degree $n$. By the universal property of the Hilbert scheme $\text{Hilb}^n_{\mathbb{P}^1_k/k} \cong \mathbb{P}^n_k$, this defines a "branch morphism", $$\text{br}_{n,\underline{\lambda}}:\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{n,\underline{\lambda}} \to \mathbb{P}^n_k.$$
<B>Proposition.</B> When the characteristic is $0$ or $p>d$ (or more generally, if the profile $\underline{\lambda}$ is tame for $p$), the branch morphism $\text{br}_{n,\underline{\lambda}}$ is étale.
<B>Proof.</B> This can be checked after base change. Thus, assume that $k$ is algebraically closed. Let $$u:C\to \mathbb{P}^1_k,$$ be a given degree $d$ $k$-morphism with $C$ a smooth, connected, projective $k$-curve of genus $g$. Let the branch points in $\mathbb{P}^1_k$ be the pairwise distinct $k$-points $(b_1,\dots,b_n)$ of $\mathbb{P}^1_k$. For each branch point, let the ramification points above $b_i$ be pairwise distinct $k$-points $(r_{i,1},\dots,r_{i,\ell_i})$ of $C$. For each $j=1,\dots,\ell_i$, let the length of $u^{-1}(b_i)$ at $r_{i,j}$ equal $\lambda_{i,j}.$ By hypothesis, each $\lambda_{i,j}$ is an integer $\geq 2$ that is prime to $p$. Thus, $\Omega_u$ is a torsion $\mathcal{O}_C$-module supported on the points $r_{i,j}$. For each $(i,j)$, denote by $\Omega_{i,j}$ the stalk of $\Omega_u$ at $r_{i,j}$. Denote by $\mathfrak{m}$ the maximal ideal of $r_{i,j}$. Then $\Omega_{i,j}$ is cyclic of length $\lambda_{i,j}-1.$ In particular, $\mathfrak{m}^{\lambda_{i,j}-1}\Omega_{i,j}$ is zero, but $\mathfrak{m}^{\lambda_{i,j}-2}\Omega_{i,j}$ is nonzero.
By the deformation theory as in the following reference,
MR0491680 (58 #10886a) <br>
Illusie, Luc <br>
Complexe cotangent et déformations. I. <br>
Lecture Notes in Mathematics, Vol. 239. <br>
Springer-Verlag, Berlin-New York, 1971. xv+355 pp. <br>
the $k$-vector space of first-order deformations of $u$ as a cover of $\mathbb{P}^1_k$ is canonically isomorphic to the $k$-vector space $H^0(C,\Omega_u)$. This, in turn, is the direct sum over all $1\leq i\leq n$ and all $1\leq j \leq \ell_i$ of the stalk $\Omega_{i,j}$. Inside this $b$-dimensional $k$-vector space, the Zariski tangent space to the component of $\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{\rho}$ is the direct sum of the submodules $\mathfrak{m}^{\lambda_{i,j}-2}\Omega_{i,j}$. This has dimension $\sum_i \ell_i$ as a $k$-vector space.
By relative duality for finite flat maps, there is a trace map of locally free $\mathcal{O}_{\mathbb{P}^1}$-modules, $$\text{Tr}_u:u_*\Omega_{C/k} \cong \textit{Hom}_{\mathcal{O}_{\mathbb{P}^1}}(u_*\mathcal{O}_C,\Omega_{\mathbb{P}^1_k/k}) \to \textit{Hom}_{\mathcal{O}_{\mathbb{P}^1}}(\mathcal{O}_{\mathbb{P}^1_k},\Omega_{\mathbb{P}^1_k/k}) = \Omega_{\mathbb{P}^1_k/k}.$$ The restriction of the trace map gives a $k$-linear isomorphism, $$\text{Tr}_{u,(i,j)}:\mathfrak{m}^{\lambda_{i,j}-2}\Omega_{i,j} \to \Omega_{\mathbb{P}^1_k/k,b_i}/\mathfrak{n}\Omega_{\mathbf{P}^1_k/k,b_i}.$$ Here $\mathfrak{n}$ is the maximal ideal at $b_i$. These $k$-linear isomorphisms define a direct sum isomorphism, $$\bigoplus_{j=1}^{\ell_i} \mathfrak{m}^{\lambda_{i,j}-2}\Omega_{i,j} \to \bigoplus_{j=1}^{\ell_i} \Omega_{\mathbb{P}^1_k/k,b_i}/\mathfrak{n}\Omega_{\mathbf{P}^1_k/k,b_i}.$$ The inverse image of the diagonal copy in the target is a $1$-dimensional $k$-subspace, $\Omega_i$. The Zariski tangent space to $\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{n,\underline{\lambda}}$ at $[u]$ is the direct sum over $i=1,\dots,n$ of the $1$-dimensional $k$-subspace $\Omega_i$. This is an $n$-dimensional $k$-vector space. Thus, $\mathcal{M}_{g,0}(\mathbb{P}^1_k/k,d)_{n,\underline{\lambda}}$ is smooth of dimension $n$. Moreover, the derivative of the branch map is the $k$-isomorphism induced by the trace maps, $$\bigoplus_{i=1}^n \Omega_i \to \bigoplus_{i=1}^n \Omega_{\mathbb{P}^1_k/k,b_i}/\mathfrak{n}\Omega_{\mathbf{P}^1_k/k,b_i}.$$
Thus, the branch map is étale. <B>QED</B>