Consider two adjacent terms in the product, and use $x=x_j,y=x_{j+1}$ for the notation. The nondecreasing inequality then gives the constraint $x \le y$ and there is a positive integer $a=j^2$ so that the contribution to the product from the $j$ and $j+1$ terms is
$$(1+a^2x^2)(1+(a+1)^2y^2). \tag{1}$$
Now if initially $x+y=k$ we only decrease (or keep the same) this two term product by replacing $x,y$ each by $k/2$, and note that since initially $x \le y$ this replacement has increased (or kept the same) $x$ while decreasing (or keeping the same) $y$, and so the replacement is consistent with the nondecreasing inequality among the $x_j$'s. It is also clear that, if $x\neq y$ were true initially, the replaced contribution would be strictly less than the initial product $(1).$
This argument shows that (as Greg Martin suggested) the minimal value of the expression occurs when all the $x_j$ are $1/n.$ If $p(n)=(2n^2+9n+1)/(6n)$ and $q(n)$ is the value of the product when each $x_j=1/n,$ then $p(1)=q(1)=2,$ so that for general $n$ the minimal lower bound of $p(n)$ happens to be exact at $n=1$. For larger $n$ it is less than $q(n)$, and has to be checked by hand for $n \le 4$, for which $q(2)-p(2)=1/4,$ and $q(3)-p(3)=53/81,$ and finally $q(4)-p(4)=655/512.$ [There may be a more clever way to see in general that $q(n)>p(n)$ but I don't see it.]
As soon as $n \ge 5$ we can use the fact that $1/n$ times the logarithm of the product is a right-endpoint sum for the integral
$$\int_0^1 \log(1+x^2)\ dx = c= \log 2+(\pi-4)/2.$$
Then from $(1/n)\cdot \log q(n) \ge c$ we get $q(n) \ge e^{cn}.$ Though I didn't do this last part carefully, it seems that for $n \ge 5$ we do have $e^{cn} >p(n)$ as required. The left side here should certainly overtake the right since it's exponential with $e^c \approx 1.302.$ I guess the reason the integral bound isn't tight for $n<5$ is because too much error occurs in replacing the sum by the integral.
EDIT (at suggestion by the OP BarbuDorel in a comment). Once it is agreed the minimum of the product occurs when all the $x_j$ are $1/n$ the product becomes
$$(1+(1/n)^2)(1+(2/n)^2)\cdots (1+(n/n)^2).$$
This is clearly at least $1$ plus the sum of the squared terms which occur as second terms in the factors. But
$$1+(1/n)^2+(2/n)^2+\cdots+(n/n)^2=\frac{2n^2+9n+1}{6n},$$
on applying the formula for the sum of the first $n$ squares. So there is no need to consider integrals as I did above. [There still remains some work to justify replacing the product of two adjacent factors $(1)$ by the product obtained on replacing adjacent $x_j$ by their average. I'll add that if I can get a simple argument.
A simpler approach. Note that the given product is bounded below by the sum
$$S=1+x_1^2+2^2x_2+\cdots + n^2x_n^2.$$ As noted above, if we can show this sum is minimal when all the $x_j$ are equal, the inequality follows. The argument for this goes similarly to before, looking at two adjacent terms, but the math is easier. Suppose the adjacent terms contribute $ax^2+by^2$ [where $0<a<b$] to the sum, and that initially $x<y$. Then if initially $x+y=k$ we have $x<k/2$. Now consider
$$ax^2+b(k-x)^2=(a+b)x^2 -2bkx+bk^2.$$
The graph of this is a parabola opening upward, with its vertex at $x=bk/(a+b)>k/2,$ where the latter follows from $0<a<b.$ So since our initial $x$ satisfies $x<k/2$ we see that the value of $ax^2+by^2$ is strictly decreased when we replace each of $x,y$ by their average $(x+y)/2.$
Now since it is clear that the restrictions on the $x_k$ define a compact bounded region in $n$ space, the function $S$ must have a minimum, and the above argument shows this minimum cannot occur at any point at which any two of the $x_j$ are not equal. Conclusion: at the minimal $S$ all the $x_j$ are equal, and the overall inequality follows from that as just outlined.