$\let\ET\bigwedge$I’ll summarize basic facts about the first question already mentioned in the comments, and add some bounds.
First, if we can write $\phi$ as $\ET_{i<j}\phi_{i,j}(x_i,x_j)$, can can expand each $\phi_{i,j}$ as a conjunction of 2-clauses (i.e., disjunctions of two literals $A,B$, which are variables or negated variables), hence boolean functions that can be written in this way are exactly those expressible by 2-CNF, aka Krom formulas. This class of formulas has attracted a lot of attention due to the fact that its satisfiability problem is efficiently solvable (in fact, NL-complete).
In order to identify formulas that express the same function, note that every Krom function $\phi\colon\{0,1\}^n\to\{0,1\}$ can be uniquely written as the conjunction of the set $C_\phi$ of all 2-clauses implied by $\phi(x_1,\dots,x_n)$. This set has the following closure properties:
(resolution) if $C_\phi$ contains $A\lor x_i$ and $B\lor\neg x_i$ for some literals $A,B$, then it also contains $A\lor B$;
(weakening) if $C_\phi$ contains $A\lor A$ (that is, $A$) for some literal $A$, it also contains $A\lor B$ for every literal $B$;
(axioms) $C_\phi$ contains all the clauses $x_i\lor\neg x_i$.
(We consider 2-clauses as unordered: $A\lor B$ is the same clause as $B\lor A$.)
Conversely, every set $C$ of 2-clauses with the three closure properties above is of the form $C_\phi$ for some Krom function $\phi$. This follows from the implicational completeness of the resolution proof system with weakening and axioms (which is even true for arbitrary clauses, not just 2-clauses).
A 2-CNF can be represented by its implication graph: a directed graph $(V,E)$ whose vertices are the literals $\{x_1,\dots,x_n,\neg x_1,\dots,\neg x_n\}$, and with a pair of edges $\neg A\to B$ and $\neg B\to A$ for every 2-clause $A\lor B$ from the 2-CNF. A graph is an implication graph of a 2-CNF iff it is skew-symmetric: $A\to B$ is an edge iff $\neg B\to\neg A$ is an edge.
The closure conditions above translate to the following conditions on the implication graph $G=(V,E)$:
$G$ is transitive;
if there is an edge from a literal $A$ to its negation, there are edges from $A$ to every vertex;
all self-loops are in the graph.
Thus, Krom functions are in 1–1 correspondence with skew-symmetric graphs with these properties. Note that conditions 1 and 3 together say that the edge relation is a preorder.
For a different characterization, we can identify a boolean function $\phi\colon\{0,1\}^n\to\{0,1\}$ with the boolean relation $\phi^{-1}(1)\subseteq\{0,1\}$. Then $\phi$ is a Krom function if and only if it has the ternary majority function
$$M(x,y,z)=(x\land y)\lor(x\land z)\lor(y\land z)$$
as a polymorphism: that is,
$$\phi(a^i_1,\dots,a^i_n)=1\text{ for $i=1,2,3$}\implies
\phi(M(a^1_1,a^2_1,a^3_1),\dots,M(a^1_n,a^2_n,a^3_n))=1.$$
As for enumeration, it is clear from the description by implication graphs that $|\mathcal P_n|$ (tabulated in OEIS: A109457 ) is related to the number of partial orders (or rather, preorders) on $n$ elements, which is
$$2^{n^2/4+O(n)}$$
by a result of Kleitman and Rothschild (their actual bound is more precise; see here for an overview).
For a lower bound, for any set of pairs $E\subseteq\binom{[n]}2$, we can form the 2-CNF
$$\phi_E(x_1,\dots,x_n)=\ET_{\{i,j\}\in E}(x_i\lor x_j)\land\ET_{i\notin\bigcup E}x_i,$$
and it is easy to see that these formulas define pairwise distinct functions. The purpose of the last conjunct is to make the function depend on all variables; that way, each $\phi_E$ gives rise to $2^n$ different functions by negating some of the variables. Thus,
$$|\mathcal P_n|\ge2^{\binom n2+n}=2^{n^2/2+n/2}.$$
For an upper bound, we need to count the number of skew-symmetric preorders $\le$ on $V=\{x_1,\dots,x_n,\neg x_1,\dots,\neg x_n\}$. I’m ignoring here condition 2, and will only use its consequence that with one exception (${\le}=V\times V$), it is not possible to have $x_i\le\neg x_i\le x_i$ for some $i$.
The standard proof that every finite partial order extends to a linear order can be easily adapted to show that every skew-symmetric partial order extends to a skew-symmetric linear order. For preorders $\le$ as above, this yields that there is a skew-symmetric linear order $\preceq$ such that
the equivalence classes of ${\approx}={\le}\cap{\le}^{-1}$ are $\preceq$-convex,
the strict order ${\le}\smallsetminus{\le}^{-1}$ is included in $\prec$.
Up to an extra factor of $2^n$, we can assume that $x_i\preceq\neg x_i$ for every $i$, thus $\preceq$ is given by a linear order on $V_0=\{x_1,\dots,x_n\}$, sitting below its reversed copy on $V_1=\{\neg x_1,\dots,\neg x_n\}$. Then, $\le$ is determined by
a preorder ${\le}_0={\le}\restriction V_0$ related to $\preceq$ as above, and
a skew-symmetric bipartite graph ${\le_1}={\le}\cap(V_0\times V_1)$, which can be identified with an undirected graph (possibly with self-loops) on $V_0$.
There are $2^n-1$ equivalence relations with $\preceq$-convex classes, and $2^{\binom n2}$ subrelations of ${\prec}\restriction V_0$, hence we obtain an elementary bound
$$|\mathcal P_n|\le2^nn!(2^n-1)2^{\binom n2}2^{\binom{n+1}2}=2^{n^2+O(n\log n)}.$$
(One could save a bit on the $2^n$ factors, but there is not much point.) If we count $\le_0$ better using the Kleitman–Rothschild bound, we obtain
$$|\mathcal P_n|\le2^{\frac34n^2+O(n)}.$$
These estimates are wasteful as a they ignore the interaction between $\le_0$ and $\le_1$, i.e., ${\le}_0\circ{\le}_1\subseteq{\le}_1$.
I have every reason to believe the correct exponent is $2^{n^2/2+\dots}$, matching the lower bound. (For example, such a bound holds when the transitive reduction of $\le_0$ is bipartite, which is true for almost all partial orders as also shown by K&R.) However, proving this seems to require an adaptation of the Kleitman&Rothschild proof, which involves an unsightly case analysis.
