I would like to answer my own question. I thank Igor Khavkine for his answer and other people for their comments which helped me to clarify the question and isolate the main difficulty.
The answer to the above question is that in Euclidean space-time the non-time-ordered (Wightman) two point function $D(\vec x,\tau):=\langle 0|\phi(\vec x,\tau)\phi(\vec 0,0)|0\rangle$ is not well defined (badly diverges) for $\tau<0$. For $\tau>0$ it is well defined and of course coincides in that domain with the time-ordered two point function. (The same is true for any derivatives of $\phi$.)
Below I will first illustrate this claim by a non-rigorous very general argument which works for not necessarily free fields, and then give explicit computation for free scalar fields.
(1) Let $\phi$ be a field, not necessarily free or scalar. The non-time-ordered product is
$$\langle 0|\phi(\vec x,\tau)\phi(\vec 0,0)|0\rangle=\langle 0|\phi(\vec x,0)e^{-\tau H}\phi(\vec 0,0)|0\rangle.$$
The spectrum of the Hamiltonian $H$ is bounded from below and unbounded from above. Generically the vector $\phi(\vec x,0)|0\rangle$
contains states of arbitrarily high energy. Thus the vector $e^{-\tau H}\phi(\vec 0,0)|0\rangle$ becomes of very infinite norm for $\tau<0$, and the two point function should badly diverge.
(2) Let now $\phi$ be a free scalar field of mass $m\geq 0$ in Euclidean space-time of dimension $D=d+1$. We exclude from the discussion the case when simultaneously $D=2$ and $m=0$. While the latter case is most important for CFT technically it requires separate treatment and will be briefly discussed below in (5). We denote by $(\vec x,\tau)$ coordinates on Euclidean space and $(\vec x,t)$ coordinates on the corresponding Minkowski space. The Wick rotation is $t=-i\tau$. On the Minkowski space-time the field $\phi(\vec x,t)$ is quantized as follows (see e.g. Peskin-Schroeder book "An introduction to QFT", Ch.2):
\begin{eqnarray}\label{1}
\phi(\vec x,t)=\int\frac{d^d\vec p}{(2\pi)^d}\frac{1}{\sqrt{2E_p}}\left(a_pe^{-ip\cdot x}+a_p^\dagger e^{ip\cdot x}\right)\big|_{p_0=E_p},
\end{eqnarray}
where $E_p=\sqrt{|\vec p|^2+m^2}$ and
\begin{eqnarray*}
[a_p,a_q]=[a_p^\dagger,a_q^\dagger]=0,\, [a_p,a_q^\dagger]=(2\pi)^d\delta^{(d)}(\vec p-\vec q),\\
a_p|0\rangle=0.
\end{eqnarray*}
The Euclidean field $\phi(\vec x,\tau)$ is obtained by substituting $t=-i\tau$. Making that and using the above commutation relations we obtain
\begin{eqnarray}\label{E:DE}
D(\vec x,\tau)=\int\frac{d^dp}{(2\pi)^d}\frac{1}{2E_p}e^{-E_p\tau+i\vec p\cdot \vec x}.(*)
\end{eqnarray}
It is obvious that for $Re(\tau)>0$ the integral absolutely converges while for $Re(\tau)<0$ it badly diverges. Notice that for purely imaginary $\tau$ it still makes sense and is the non-time-ordered two point function in the Minkowski space-time.
(3) One can show that the function $D(\vec x,\tau)$ in the formula ($*$) in the half-space $\tau>0$ is $SO(D)$-invariant, i.e. depends only on $|\vec x|^2+\tau^2$ (e.g. for $D=2$ it is shown in Rotation invariance of an integral ).
Let us show that this implies that the scalar free field does not admit a quantization which is equivariant under the whole Poincare group, i.e. all orientation preserving isometries of $\mathbb{R}^D$. This is in sharp contrast with the Minkowski case, and I guess this is what was meant by Igor Khavkine when he said that there is no operator formalism in Euclidean case. Indeed otherwise for any given point in $x\in\mathbb{R}^D\backslash\{0\}$ we always could choose a time direction such that such that the time coordinate would be positive: $x=(\vec x,\tau), \tau>0$. Then $\langle 0|\phi(x)\phi(0)|0\rangle$ would be well convergent, this is a contradiction. However my impression is that if one fixes the time direction then apparently the free scalar field can be quantized in a way equivariant under the smaller group of all translations and rotations preserving this fixed time direction.
(4) The formula ($*$) defines the function $D$ which is $SO(D)$-invariant in the half-space $\tau>0$ and clearly satisfies $(\Delta-m^2)D=0$. By rotation invariance we can extend $D$ uniquely to $\mathbb{R}^D\backslash\{0\}$. This extension satisfies the same equation. One may ask what is the physical meaning of this extended function? Obviously it is the time-ordered two point function which well converges everywhere but $0$.
(5) Let us discuss briefly the case $D=2,m=0$. It is clear that the integral ($*$) has a divergency at $p=0$ in this case (as opposed to all previously discussed cases). Already in the Minkowski 2d space-time there is an (infrared) divergency which makes the above procedure of quantization impossible. It was shown by Wightman ("Introduction to some aspects of quantizes fields", in "Lectures notes, Cargese Summer School, 1964") that one has to abandon some of the usual requirements imposed on quantization; he has shown that it is possible to quantize the free massless scalar field in 2d Minkowski space-time provided one abandons the assumption of positivity of the scalar product in Hilbert space. I did not try to check how to apply the Wick rotation as in (2) to Wightman's procedure of quantization.
In CFT one works with the periodic boundary conditions: $\phi(x+L,\tau)=\phi(x,\tau)$. Isometries (and even conformal transformations) of the cylinder are just translations preserving the time direction. Thus the possibility of operator formalism of quantization does not contradict to what was said by Khavkine due to the discussion in (3).
