It seems that the output of sklearn.metrics.pairwise.euclidean_distances is different to the formula on doc
The sklearn docs' formula says it is applying to row vectors $x$ and $y$.  When you call np.dot on the matrices $X$ and $Y$ it takes the matrix product.

EDIT (responding to question in comments):
It's not straightforward, as the row-vs-row operations needed aren't quite the usual matrix operations.  The source code for euclidean_distances does it this way (except that it does lots of input checks, operates on sparse inputs when possible, etc.):

(X*X).sum(axis=1)[:, np.newaxis] - 2*np.dot(X,Y.T) + (Y*Y).sum(axis=1)[np.newaxis, :]


That's not exactly straightforward itself, so I'll say a little more.  Say $X$ has $m$ rows and $Y$ has $n$ rows.  The middle term, by taking $Y^T$, gives us a $m\times n$ matrix whose $(i,j)$-entry is the dot product of the $i$th row of $X$ with the $j$th row of $Y$.  In the other terms, * on numpy arrays is the coordinate-wise product; summing along rows gives us the rows' squared-norms.  The newaxis is a nice trick: casting the first term to now be a $m\times 1$ matrix, adding it to the middle term's $m\times n$ matrix actually adds it to every column of that matrix (without needing to actually build the matrix of repeated columns of $X$'s squared-norms).  And of course similarly for the last term: casting to a $1\times n$ matrix makes it add to every row of the result.