Let $\theta \in \Theta \subset \mathbb R^m$ and $\phi : \Theta \to \mathbb R^d$ so $T(x) \in \mathbb R^d$ too. I'll use without proof the fact that differentiation and integration can be exchanged for exponential families.
First I'm going to treat this as a natural exponential family and see what happens (this is like if $\phi$ is just the identity function). I'll use $\theta$ as the parameter so right now I have $p(x|\theta) = \exp(\theta^TT(x) - A(\theta))f(x)$.
In this case the differentiation is easy and we have
$$
0 = \int p(x|\theta) \left[T_j(x) - \nabla _jA(\theta)\right]\,\text dx
$$
so
$$
\nabla A(\theta) = \text E(T(X)).
$$
So if you're content to work with $\phi := \phi(\theta)$ as the parameter, then this is a tidy result.
Now I'll consider what happens if $m = d = 1$ to build our intuition for the general case. We'll have
$$
0 = \int p(x|\theta)\left[T(x) \cdot \phi'(\theta) - A'(\theta) \right]\,\text dx \\
\implies \text E(T(X)) \cdot \phi'(\theta) = A'(\theta)
$$
so
$$
\text E(T(X)) = \frac{A'(\theta)}{\phi'(\theta)}
$$
as you thought.
But the situation is more complicated if $m,d > 1$. The derivates in particular are more complicated because now $\phi$ is vector-valued so we'll get a matrix of first derivatives.
We need
$$
\frac{\partial}{\partial \theta_j} \phi(\theta)^TT(x).
$$
Writing this as a sum, we have
$$
\frac{\partial}{\partial \theta_j} \phi(\theta)^TT(x) = \sum_{i=1}^d T_i(x) \frac{\partial}{\partial \theta_j}\phi_i(\theta).
$$
$\phi : \Theta\to\mathbb R^d$ with $\Theta\subset\mathbb R^m$ so for each component we have $\phi_i : \Theta\to\mathbb R$ with $i=1,\dots,d$. Each $\phi_i$ is just a typical scalar-valued function so $\frac{\partial}{\partial \theta_j} \phi_i = (\nabla \phi_i)_j$.
I'll collect these gradients into a matrix $\Phi(\theta) \in \mathbb R^{m\times d}$ where
$$
\Phi(\theta) = \left[\begin{array}{c|c|c|c}\nabla \phi_1(\theta) & \nabla \phi_2(\theta) & \cdots & \nabla \phi_d(\theta)\end{array} \right].
$$
This means
$$
\frac{\partial}{\partial \theta_j}\phi_i(\theta) = \Phi(\theta)_{ji}
$$
so
$$
\frac{\partial}{\partial \theta_j} \phi(\theta)^TT(x) = \left(\Phi(\theta)T(x)\right)_j.
$$
All together this means
$$
\mathbf 0 = \int p(x|\theta) \left[ \Phi(\theta)T(x) - \nabla A(\theta) \right] \,\text dx \\
\implies \Phi(\theta)\text E [T(X)] = \nabla A(\theta)
$$
which is a linear system that we're trying to solve for $\text E [T(X)]$. The dimensions in question and the specific properties of $\phi$ now will determine how many solutions this system has.
Here's one example. Suppose $X \sim \mathcal N(\mu, \sigma^2)$ so
$$
p(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).
$$
Rearranging and expanding the quadratic we have
$$
p(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi}}\exp\left(\frac{\mu}{\sigma^2}x - \frac{1}{2\sigma^2} x^2 - \left[\frac{\mu^2}{2\sigma^2} + \frac 12 \log \sigma^2\right]\right)
$$
which means we'll have $\theta = (\mu,\sigma^2)$, $\Theta = \mathbb R \times (0,\infty)$, $T(x) = {x \choose x^2}$,
$$
\phi(\theta) = \left(\begin{array}{c}\frac{\theta_1}{\theta_2} \\ \frac{-1}{2\theta_2}\end{array}\right),
$$
and
$$
A(\theta) = \frac 12 \theta_1^2\theta_2^{-1} + \frac 12 \log\theta_2.
$$
This means
$$
\Phi(\theta)= \left(\begin{array}{cc}\theta_2^{-1} & 0 \\ -\theta_1\theta_2^{-2} & \frac 12 \theta_2^{-2}\end{array}\right).
$$
Noting that
$$
\det \Phi(\theta) = \frac 12 \theta_2^{-3} > 0
$$
we know this matrix is invertible, and it's just $2\times 2$ so the inverse is easy enough to compute:
$$
\Phi^{-1}(\theta) = 2\theta_2^3 \left(\begin{array}{cc}\frac 12 \theta_2^{-2} & 0 \\ \theta_1\theta_2^{-2} & \theta_2^{-1}\end{array}\right) = \left(\begin{array}{cc}\theta_2 & 0 \\ 2\theta_1\theta_2 & 2\theta_2^2\end{array}\right).
$$
We also have
$$
\nabla A(\theta) = \left(\begin{array}{cc}\theta_1\theta_2^{-1} \\ -\frac 12\theta_1^2\theta_2^{-2} + \frac 12 \theta_2^{-1}\end{array}\right)
$$
so
$$
\text E(T(X)) = \Phi^{-1}(\theta)\nabla A(\theta) \\
= \left(\begin{array}{c}\theta_1 \\ \theta_1^2 + \theta_2\end{array}\right) = {\mu \choose \mu^2 + \sigma^2}
$$
which we know to be the first two raw moments of $X$. This was a lot of work to get there but it shows the unavoidable role of the matrix $\Phi$.
Update: here's the multivariate normal case. Let $X \sim \mathcal N_n(\mu, \Omega^{-1})$ where $\Omega$ is the precision matrix. I'll use $\Omega_{(i)}$ for the $i$th row of $\Omega$, and for a vector $v$ and $1\leq i\leq j$ I'll use $v_{i:j}$ to be $(v_i,v_{i+1},\dots,v_j)$.
The density is
$$
p(x \mid \mu, \Omega) = \exp\left[-\frac 12 x^T\Omega x + x^T \Omega \mu - \frac 12 \left(\mu^T\Omega\mu - \log\det \Omega\right)\right] \cdot(2\pi)^{-n/2}.
$$
The parameters are $\mu$ and $\Omega$ but I want to deal with these as a flat vector so I'll combine $\mu$ with $\Omega$ "unrolled" as
$$
\theta = (\mu, \Omega_{11}, \Omega_{12}, \dots, \Omega_{1n}, \Omega_{22}, \dots, \dots, \Omega_{n-1,n}, \Omega_{nn}).
$$
I'm only taking the diagonal and one triangle of $\Omega$ since it's symmetric. This means $\Theta\subset \mathbb R^m$ where $m = n + {n+1\choose 2} = 2n + {n\choose 2}$.
I now need to get this as $\phi(\theta)^TT(x) - A(\theta)$. For $\phi$ and $T$ I'll use the fact that
$$
x^T\Omega x = 2\sum_{i < j} x_ix_j \Omega_{ij} + \sum_{i=1}^n x_i^2\Omega_{ii}
$$
so I'll take
$$
\phi(\theta) = \left(\Omega\mu, -\frac 12 \Omega_{11}, -\Omega_{12}, \dots, -\Omega_{1n}, -\frac 12 \Omega_{22}, -\Omega_{23}, \dots, \dots, -\Omega_{n-1,n}, -\frac 12 \Omega_{nn}\right)
$$
(I'm writing this in terms of $\mu$ and $\Omega$ instead of $\theta$ just because it's way clearer, but these can all be written in terms of $\theta$ if desired) and
$$
T(x) = (x, x_1^2, x_1x_2, \dots, x_1x_n, x_2^2, x_2x_3, \dots,\dots, x_n^2).
$$
We can now construct $\Phi$. I'll build it as a $2\times 2$ block matrix where
$$
\Phi(\theta) = \left(\begin{array}{c|c}A&B\\\hline C&D\end{array}\right)
$$
with $A$ being $n\times n$ and $D$ being ${n+1c\choose 2}\times{n+1\choose 2}$. $A$, $B$, and $D$ are easy to work out.
For $A$, I need $\frac{\partial \phi_i}{\partial \theta_j}$ where $i,j=1\dots,n$. Since $\phi_i(\theta) = \Omega_{(i)}^T\mu$ for this range of $i$ values, and $\theta_j = \mu_j$, we'll get
$$
(\nabla \phi_i)_{1:n} = \Omega_{(i)}
$$
(using the symmetry of $\Omega$ for the fact that $\Omega_{k\ell} = \Omega_{\ell k}$). This means that $A = \Omega$.
Now for $i > n$, $\phi_i(\theta)$ is proportional to some $\Omega_{j\ell}$ so derivatives w.r.t. $\mu$ are zero. This means $B = \mathbf O$, the appropriately sized zero matrix.
Next, for $i,j>n$ if $i\neq j$ then $\frac{\partial \phi_i}{\partial \theta_j} = 0$ so $D$ is diagonal. If $i$ corresponds to a diagonal entry of $\Omega$ then we have a derivative of $-\frac 12$, otherwise it'll be a derivative of $-1$. The first element of the diagonal of $D$ corresponds to $\Omega_{11}$ so it begins with $-\frac 12$. The next value of $-\frac 12$ is for $\Omega_{22}$ which happens at element $1+n$ since all of the first row of $\Omega$ has been counted and $\Omega_{21}$ is skipped (since it's accounted for by $\Omega_{12}$). There are $n-2$ remaining elements of $\Omega_{(2)}$ to account for, so the 3rd value of $-\frac 12$ in $D$ is at position $1+n+(n-1)$. The next $-\frac 12$ is at $1+n+(n-1)+(n-2)$ and so on, so the locations of $-\frac12$s are given by the sequence
$$
1+\sum_{j=0}^{i-2}(n-j)
$$
with $i=1,\dots,n$. As an example, if $n=5$ then the $-\frac 12$s will be at locations $1,6,10,13,15$. As a sanity check, $D$ is supposed to have a diagonal of length ${n+1\choose 2}$, and the final $-\frac 12$ is for $\Omega_{nn}$. We have
$$
1+\sum_{j=0}^{n-2}(n-j) = \sum_{j=0}^{n-1}(n-j) \\
= \sum_{j=1}^n j = {n+1\choose 2}
$$
so this agrees with that.
Finally, we need $C$. I'm going to skip $C$ for now until the final step.
For $\Phi^{-1}$, $\Omega$ and $D$ are always invertible so $\Phi$ is too, and we can use the formula for the inverse of a $2\times 2$ block matrix to find
$$
\Phi^{-1}(\theta) = \left(\begin{array}{c|c}\Omega^{-1} & \mathbf O \\ \hline
-D^{-1}C\Omega^{-1} & D^{-1}\end{array}\right).
$$
Now we can get $\nabla A$ to finish this off. I'll also partition it into the first $n$ elements, which lines up the upper two blocks of $\Phi^{-1}$, and the remaining ${n+1\choose 2}$ which I'll call $V$. We have
$$
A(\theta) = \frac 12 \mu^T\Omega\mu - \frac 12 \log\det \Omega
$$
so
$$
\frac{\partial A}{\partial \theta_{1:n}} = \frac{\partial A}{\partial \mu} = \Omega\mu.
$$
This means
$$
\nabla A(\theta) = {\Omega\mu \choose V}
$$
so
$$
\Phi^{-1}(\theta)\nabla A(\theta) = {\mu \choose -D^{-1}(C\mu - V)}.
$$
This shows that $\text E(T_{1:n}(X)) = \text E(X) = \mu$ which is nice to see.
For the rest, I won't work it out in general but I'll just show what we get for the elements corresponding to diagonal terms of $\Omega$ and off-diagonal terms.
For $V$, I'll need derivatives of $\log\det \Omega$ w.r.t. $\Omega_{ij}$ which I can get via this result:
$$
\frac{\partial}{\partial \Omega_{ij}} \log\det \Omega = (2 \Sigma - \text{diag}(\Sigma))
$$
due to symmetry. This means that for the elements of $V$ corresponding to a diagonal element of $\Omega$ we'll have $\frac 12 \mu_i^2 - \frac 12 \Sigma_{ii}$ while the off-diagonal will be $\mu_i\mu_j - \Sigma_{ij}$.
For $C\mu$, the diagonal elements contribute $\mu_i^2$ because for that row the only non-zero element is $\mu_i$ in the $i$th column. $-D^{-1}$ will multiply the result by $2$, and then $V_i$ is
$$
2\left(\mu_i^2 - \frac 12 \mu_i^2 + \frac 12 \Sigma_{ii} \right) \\
= \mu_i^2 + \Sigma_{ii}.
$$
For the off-diagonal, you can show that $C\mu$ will contribute $2\mu_i\mu_j$ so we get
$$
2\mu_i\mu_j - \mu_i\mu_j + \Sigma_{ij} = \mu_i\mu_j + \Sigma_{ij}.
$$
Thus after all this work we have confirmed that $\text E(X_i) = \mu_i$ and $\text E(X_iX_j) = \mu_i\mu_j + \Sigma_{ij}$.
