Let $E$ denote the interior of the ellipse and $\partial E$ the boundary (i.e. the ellipse itself). You're trying to calculate
$$\color{blue}{\oint_{\partial E}\mathbf{F}\cdot d\mathbf{r}}$$
Your first hope is that Green's theorem applies, so you can instead calculate
$$\int_E(Q_x-P_y)\,dA$$
The problem is that the hypotheses of Green's theorem are not satisfied, for the functions $P$ and $Q$ are not continuously differentiable throughout the region $E$. As you observed, there is no way we may continuously extend the field to be defined at the point $(1,0)$, so that point is a singularity.
The way to proceed in such cases is to put a little closed curve around the singularity but still lying entirely inside the region $E$ of interest; in this case a unit circle will be convenient (the calculation will show you why). So let $D$ be the unit disk centered at $(1,0)$, and let $\partial D$ be the unit circle centered there, which has parametrization $C(t)=(\cos t +1, \sin t)$ for $t\in[0,2\pi)$.
You can easily check by direct calculation that
$$\color{red}{\oint_{\partial D}\mathbf{F}\cdot d\mathbf{r}}=\int_0^{2\pi}(\cos^2t+\sin^2t)\,dt=2\pi$$
Now, consider the region $E-D$, the complement of the disk in the elliptical region. In this region, $P$ and $Q$ are continuously differentiable everywhere, so we may use Green's theorem, which says
$$\int_{E-D}(Q_x-P_y)\,dA=\oint_{\partial(E-D)}\mathbf{F}\cdot d\mathbf{r}\tag{$\star$}$$
The left side vanishes, because $Q_x-P_y=0$ in the region $E-D$. What about the right side? What is the "boundary" $\partial(E-D)$ of the elliptical annular region $E-D$? It is $\partial E - \partial D$, where we interpret the negative sign on $\partial D$ to mean reversing the standard counterclockwise orientation.
So $(\star)$ says
$$0=\color{blue}{\oint_{\partial E}\mathbf{F}\cdot d\mathbf{r}}-\color{red}{\oint_{\partial D}\mathbf{F}\cdot d\mathbf{R}}$$
The blue integral is the one the problem asks for, and the red integral we calculated above.
To summarize, the strategy for calculating line integrals of curl-free fields around closed curves that encircle a single singularity is to replace the original closed curve with an easier one (but still contained within the original one). Green's theorem in the region between the curves guarantees the line integral around the easier curve equals the line integral around the original curve.
What makes a good candidate for an "easier curve"? Well, you're trying to find a curve whose parametrization will make the integrand simple enough to calculate. In your problem, integrating around the ellipse gives a "hard" integrand. But integrating around a circle centered at the singularity gives an "easy" integrand.
Why a circle in this case? Because your field is closely related to the pullback of the length element of a circle of radius $r$ (the length element is $ds=r d\theta$); thus it was, in some sense, tailor-made to be integrated around a circle. It is not always so easy to spot what sort of curve will make the line integral easier, but you should remember the general form $(-\frac{y}{r^2}, \frac{x}{r^2})$, which is $\frac{1}{r}$ times the circle's unit tangent vector. The line integral of the unit tangent vector around the circle gives $2\pi r$, and the factor $\frac{1}{r}$ scales the result to $2\pi$. (In other words, you are just integrating $d\theta=\frac{1}{r}ds$ around the circle.) So the radius of the circle you use around the singularity doesn't matter. (This make sense, because the region between two concentric circles centered at the singularity doesn't contain a singularity and the field is conservative there, so by the above argument, the line integral of the field around them should be the same.)
Any mapping, be it a vector field or a scalar function or something else, requires a domain.
It is true that where $\phi(x,y)$ is defined, $\nabla \phi = \vec F$. But $\vec F$'s domain is the plane minus the origin, while $\phi$'s domain is the plane minus a line (the $y$-axis).
Since there's no function with the same domain as $\vec F$ whose gradient is $\vec F$, $\vec F$ is not conservative.
Notice that the right half of the plane is simply connected, and as you've shown, $\vec F$ restricted to that domain is conservative. $\phi$ works as a potential on that domain.
The upshot is that the question of whether $\vec F$ is conservative on $U$ is a question not just about the component functions of $\vec F$ but the “shape” (we say topology) of $U$.
Best Answer
One of the most important "functions" in classical analysis is the ${\rm arg}$ function, giving the polar angle of a point $(x,y)$, resp., of $z=x+iy$, in the punctured plane "up to an additive constant $2k\pi$". Locally, i.e., in suitable neighborhoods of points $(x_0,y_0)\in\dot{\mathbb R}^2$ this function has well-defined real representants. E.g., in the half plane $x>0$ we may choose the representant $$\phi(x,y):=\arctan{y\over x}\ .$$ Now this ${\rm arg}$ "function" has a well-defined gradient $$\nabla{\rm arg}(x,y)=\left({-y\over x^2+y^2},{x\over x^2+y^2}\right)\qquad\bigl((x,y)\ne(0,0)\bigr)\ .$$It turns out that your ${\bf F}$ is nothing else but $\nabla{\rm arg}$. Since, locally, ${\rm arg}$ is represented by smooth functions $(x,y)\mapsto\phi(x,y)$ it follows that ${\rm curl}\,{\bf F}\equiv0$.
But it is impossible to concatenate the local polar angles $\phi(x,y)$ to one single real-valued polar angle function $\phi:\>\dot{\mathbb R}\to{\mathbb R}$. This then implies that the given ${\bf F}$ is not conservative. A definitive proof of this fact comes from computing $\int_\gamma {\bf F}\cdot d{\bf z}$ along the unit circle $\gamma$. Along this circle the argument increases continuously by $2\pi$, hence the value of this integral is $2\pi\ne0$.