The functions which preserve asymptotic equivalence have turned out to be one of a vast number of generalizations of Karamata's regularly varying functions. Paper "On some extensions of Karamata's theory and their applications" by Buldygin et. al. introduces pseudo-regularly varying (PRV) functions which are, in fact, exactly that measurable functions, not necessarily tending to infinity, which preserve asymptotics at $+\infty$ in my sense.
The paper contains many properties of PRV functions and conditions for a function to be PRV. For example, PRV fuctions are characterized as being eventually ($\forall x>x_0$) of the form
$$
f(x) = \exp \bigg( a(x) + \int\limits_{x_0}^x b(t) \frac{dt}{t} \bigg),
$$
where $a$ and $b$ are bounded measurable functions with $\lim\limits_{c\rightarrow 1}\limsup\limits_{x\rightarrow +\infty} \lvert a(cx)-a(x) \rvert = 0$.
So the reference request is fulfilled, I'll return to finer class $\mathcal A$ described in the question and answer the "mini-questions".
Besides the above, we can charactirize $\mathcal A$ in terms of uniform continuity, although this characterization is not so constructive.
Proposition. A continuous function $f\colon \mathbb R_+ \rightarrow \mathbb R_+$ with $\lim \limits_{x \rightarrow +\infty} f(x) = +\infty$ belongs to $\mathcal A$ if and only if function $F(t) = \log f(e^t)$, correctly defined on some ray $[T,+\infty)$, is uniformly continuous.
Proof. Note that $\log f(x) = F(\log x)$.
1) To prove the "if" part, pick $\{a_n\}_{n=1}^\infty$ and $\{b_n\}_{n=1}^\infty$ as in the definition of $\mathcal A$. Since $\lim\limits_{n\rightarrow\infty} \log\tfrac{a_n}{b_n} = 0$ and $F$ is uniformly continuous,
\begin{align*}
\left| \log \frac{f(a_n)}{f(b_n)} \right| &= \lvert \log f(a_n) - \log f(b_n) \rvert = \lvert F(\log a_n)-F(\log b_n) \rvert = \\ &= \lvert F(\log b_n+\log\tfrac{a_n}{b_n}) - F(\log b_n) \rvert \leqslant \sup\limits_{t \geqslant T} \lvert F(t+\log\tfrac{a_n}{b_n}) - F(t) \rvert \xrightarrow{\;n\rightarrow\infty\;} 0,
\end{align*}
so finally $f(a_n) \sim f(b_n)$ as $n$ tends to infinity.
2) To prove the "only if" part, suppose that, coversely, $F$ is not uniformly coninuous. That means there exists $\varepsilon > 0$ and two sequences $\{t_n\}_{n=1}^\infty, \{s_n\}_{n=1}^\infty \subseteq [T, +\infty)$, such that
$$
\lim\limits_{n\rightarrow\infty} (t_n-s_n) = 0 \quad\text{and}\quad \forall n\in \mathbb N\ \; \lvert F(t_n) - F(s_n) \rvert \geqslant \varepsilon.
$$
It's necessary that $\{t_n\}_{n=1}^\infty$ is unbounded, because by Cantor's theorem $F$ is uniformly continuous on each finite interval. Hence, taking subsequence if needed, we can assume that $t_n$ (and so $s_n$) goes to infinity.
Define $a_n := \exp(t_n), \; b_n := \exp(s_n)$. We have $\lim\limits_{n\rightarrow\infty} \tfrac{a_n}{b_n} = \lim\limits_{n\rightarrow\infty} \exp(t_n-s_n) = 1$ and
$$
\bigg\lvert \! \log \frac{f(a_n)}{f(b_n)} \! \bigg\rvert = \lvert F(\log a_n)-F(\log b_n) \rvert = \lvert F(t_n)-F(s_n) \rvert \geqslant \varepsilon \;\;\; \forall n \in \mathbb N,
$$
what contradicts with $f \in \mathcal A$.
From this characterization it follows that a functions from $\mathcal A$ grows at most like power function. (It's also follows from the exponent expression for PRV functions.)
Corollary. If $f \in \mathcal A$, then $f(x) = O(x^\alpha)$ for some $\alpha > 0$ as $x \rightarrow +\infty$.
Proof. Take $F$ like in above proposition. Since $F$ is uniformly continuous, $F(t) \leqslant \alpha t \;\;\; \forall t \geqslant T$ with some $\alpha,T \geqslant 0$ (see this question). Then
$$
f(x) = \exp\big(F(\log x)\big) \leqslant \exp(\alpha \log x) \leqslant x^\alpha \quad \forall x \geqslant e^T.
$$
And finally, we can indeed replace equivalent sequences in the definition of $\mathcal A$ with such $C^\infty$-smooth functions on $\mathbb R_+$ and get the same class. It's straightforward, we just need to construct smooth function with given values in, say, $\mathbb N$.
We may also restrict the sequences from the definition to be strictly increasing, hence each real sequence not bounded above admits a strictly increasing subsequence. And thus the smooth functions from the alternative definition can be all stricly increasing. It's bit more tricky and this is crucial for l'Hopital's rule application to prove
Proposition. If $f, g \in \mathcal A$, then $F(x)=\int_0^x f(t)dt$ and $G(x) = \int_0^x g(t)dt$ lie in $\mathcal A$ and, moreover, $F \sim G$ at infinity.
Following a more elementary route we will first break the integral in two pieces:
$$B(\frac{N+1}{N+2};N+1,p+1)=B(N+1,p+1)-I\\I=\int^{1}_{\frac{N+1}{N+2}}x^{N}(1-x)^pdx$$
Perform the following change of variables $x=1-\frac{t}{N+2}$:
$$I=\frac{1}{(N+2)^{p+1}}\int_{0}^1t^p(1-\frac{t}{N+2})^{N}dt$$ However it is possible to estimate that
$$(1-\frac{t}{N+2})^{N}=e^{-t}(1+\frac{4t-t^2}{2N}+\mathcal{O}(N^{-2}))$$
and performing the integrals we get:
$$I\sim\frac{1}{N^{p+1}}\Big[\gamma(p+1,1)+\frac{2}{N}(\gamma(p+2,1)-p-1-\frac{1}{4}\gamma(p+3,1))\Big]$$
While on the other hand we may write
$$B(N+1,p+1)\sim \frac{\Gamma(p+1)}{N^{p+1}}\Big[1-\frac{(p+1)^2+3(p+1)}{2N}\Big]$$
So all in all as a leading approximation one gets:
$$B(\frac{N+1}{N+2};N+1,p+1)\sim\frac{\Gamma(p+1)-\gamma(p+1,1)}{N^{p+1}}=\frac{\Gamma(p+1,1)}{N^{p+1}}+\mathcal{O}(N^{-p-2})$$
where $\Gamma(z,x)=\int_x^{\infty}t^{z-1}e^{-t}dt$ is the upper incomplete Gamma function. From the above one can calculate also the subleading term. More subleading terms can also be readily calculated but probably not in closed form.
Best Answer
Your particular result is true because it happens to be a Laplace Transform of a function that does not grow too quickly near the postiive real line.
Indeed, your book of Gamelin (page 367) expresses the trigamma $\psi'$ in this form,
$$ \psi'(z)=\frac{d}{dz} \frac{\Gamma'(z)}{\Gamma(z)} = \int_0^\infty e^{-zt}\frac{t}{1-e^{-t}} dt, \quad \Re z> 0.$$ This works because $$f(t) = \frac{t}{1-e^{-t}}$$ is exponential type of order $0+$ on a small neighbourhood of the positive real line, which fits into a more general result:
Implicit above is the fact that $\int_0^\infty e^{-zt}f(t) dt$ converges for analytic functions of exponential order $\rho$ in $U_\delta$. We say $f$ is of exponential order $\rho>0$ in $D$ if there exists $M$ such that for every $t\in D$, $$ |f(t)| \le M e^{\rho |t|}.$$
An example: if $f(t)=1$ then we can take any $\rho>0$ and $\mathcal Lf = \frac{1}{z}\simeq \frac1z + 0 + \dots $ is indeed valid on every half plane $ \{\Re z > \rho' \}$.
The proof is relatively straightforward; here's the sketch. The appearance of derivatives strongly suggests integration by parts. Under the given assumptions, you can check that $f^{(k)}$ is also of exponential type $\rho$ for every $k$, which justifies the integration by parts. We record the constant that verifies this as $M_k$, i.e. $$ |f^{(k)}(z)| \le M_k e^{\rho |z|}.$$
Performing the integration by parts inductively leads to
$$\mathcal{L} f(z)=\sum_{k=0}^{n-1} \frac{f^{(k)}(0)}{z^{k+1}}+\frac{1}{z^n} \int_{0}^{\infty} e^{-z t} f^{(n)}(t) d t$$ So we need to bound this last integral. One last integration by parts gives \begin{align} \left| \int_{0}^{\infty} e^{-z t} f^{(n)}(t) d t\right| & \le \frac{|f^{(n)}(0)|}{|z|} + \frac1{|z|}\int_0^\infty \underbrace{|e^{-tz} f^{(n+1)} (t)|}_{\Large\le M_{n+1} e^{-t \Re z+\rho t}} dt \\ &\le \frac{M_n}{|z|} + \frac{M_{n+1}}{|z|(\Re z-\rho)} \\ &\le \frac{M_n}{|z|} + \frac{M_{n+1}}{|z|(\rho'-\rho)} \\&=: \frac{C_{n+1}}{|z|} \end{align} where $\rho<\rho'<\Re z$ was crucially used. This gives the required estimate to satisfy the definition of an asymptotic series (with asymptotic sequence $(z^{-k-1})_{k\ge 0}$).
Here is the bibliography from the description of my old course on Asymptotic Analysis, may be useful...