My question is:
How can I get the ditribution of $\frac{S_{t+1}-S_{t}}{ S_{t}}$ ? (at least its expectation and standard deviation under filtration $F_t$)

If we come back to the first equation it seams that if $\mathrm{d}t$ is small enough then it is a normal distribution with mean $k(\theta - \ln S_{t}) \mathrm{d}t$ and sd $\sigma\mathrm{d}t$, but I want to find mathematically if it is true, and if yes, under what assumptions.

1 ) A first-order Taylor expansion gives $\ln \left(\frac{S_{t+\Delta_t}}{S_t}\right)\approx \frac{S_{t+\Delta_t}-S_{t}}{S_t}+o(\Delta_t)$ , thus unless $\Delta_t$ is not small you can drop the residual term and consider $Z_t\overset{law}{=}\frac{S_{t+\Delta_t}-S_{t}}{S_t}$.

2 ) Calculation of the moments: we can proceed by using the classical Dynkin way

we remind that $dZ_t=k(\theta-\frac{\sigma^2}{2k}-Z_t)dt+\sigma dW_t$, so

so, your p-th moment is the solution of the above deterministic ODE which can be solved step-by-step starting $p=1$.

3 ) Exact distribution of the yield: in general it's far to be simple to get exact distribution/ simulation for diffusions (and more for multidimensional SDE's -- example: Heston's joint SDE has an exact pdf which computation requires Malliavin calculus).

So, in my opinion unless you are considering simulations with large time step it would be useless to deal with the exact distribution of your SDE.