Your substitution is not correct. You should still still have a $\sin^3 x$ upstairs.
However when integrating a product of an even power of $\tan$ with an even power of $\sec$, you can do the following, which takes advantage of the facts that $\tan^2x+1=\sec^2 x$ and that the derivative of $\tan x$ is $\sec^2 x$:
First write
$$\eqalign{\tan^4 x\,\sec^6 x&=\tan^4x\,\sec^4x\cdot \sec^2 x\cr
&= \tan^4x \,(\sec^2x)^2\cdot \sec^2 x\cr
&=\tan^4x \,(\tan^2 x+1)^2\cdot \sec^2 x.}$$
To evaluate the integral, set $u=\tan x$ (so that $du=\sec^2 x \,dx$).
We'll consider the integral
$$
J(x) = \int_0^\infty \frac{1}{1+t} e^{-xt^2}\,dt,
$$
which is in the form
$$
\int_0^\infty f(t) e^{xg(t)}\,dt
$$
with
$$
f(t) = \frac{1}{1+t} \qquad \text{and} \qquad g(t) = -t^2.
$$
The basic idea of the Laplace method is that the main contribution to the integral $J(x)$ comes from the points $t$ where the integrand $f(t)e^{xg(t)}$ is largest.
As $x$ grows larger, the behavior of the factor $f(t)$ matters less and less; where $g(t)$ is large $f(t)e^{xg(t)}$ will be large, and where $g(t)$ is small $f(t)e^{xg(t)}$ will be small. This means that we're really interested in the maxima of $g(t)$ -- these will be the points $t$ where the integrand $f(t)e^{xg(t)}$ is largest when $x$ is large.
For this problem the function $g(t) = -t^2$ has a single maximum at $t=0$. By the heuristics aboove, the large-$x$ behavior of $J(x)$ should be very similar to the large-$x$ behavior of an integral over a small neighborhood of this point:
$$
J(x) \approx \int_0^\delta \frac{1}{1+t} e^{-xt^2}\,dt.
$$
Our goal is to justify this and to further simplify the integrand to get a simplified expression for the approximation.
The first thing we want to do is show that the integral over the region $[\delta,\infty)$ is very small. Here we can get the estimate
$$
\begin{align}
0 \leq \int_\delta^\infty \frac{1}{1+t} e^{-xt^2}\,dt &< \frac{1}{1+\delta} \int_\delta^\infty e^{-xt^2}\,dt \\
&= \frac{1}{1+\delta} e^{-\delta^2x} \int_\delta^\infty e^{x (\delta^2-t^2)}\,dt \\
&\leq \frac{1}{1+\delta} e^{-\delta^2x} \int_\delta^\infty e^{\delta^2-t^2}\,dt \tag{1}
\end{align}
$$
for all $x \geq 1$. Thus
$$
J(x) = \int_0^\delta \frac{1}{1+t}e^{-xt^2}\,dt + O\!\left(e^{-\delta^2x}\right) \tag{2}
$$
for all $x \geq 1$.
Our next goal is to study this remaining integral and to show that, to leading order, it is larger than the error term $O(e^{-\delta^2x})$. To do this we will use the fact that
$$
\lim_{t \to 0} \frac{1}{1+t} = 1.
$$
Let
$$
h(t) = \frac{1}{1+t} - 1 = \frac{-t}{1+t},
$$
so that
$$
\int_0^\delta \frac{1}{1+t} e^{-xt^2}\,dt = \int_0^\delta e^{-xt^2}\,dt + \int_0^\delta h(t) e^{-xt^2}\,dt. \tag{3}
$$
We expect the second integral to be smaller than the first since $h(t) \approx 0$ for small $\delta$. Indeed,
$$
\begin{align}
\int_0^\delta e^{-xt^2}\,dt &= \int_0^\infty e^{-xt^2}\,dt - \int_\delta^\infty e^{-xt^2}\,dt \\
&= \frac{1}{2} \sqrt{\frac{\pi}{x}} + O\!\left(e^{-\delta^2 x}\right), \tag{4}
\end{align}
$$
where we used an argument very similar to $(1)$ to obtain the $O(\cdots)$ estimate, and
$$
\begin{align}
\left| \int_0^\delta h(t) e^{-xt^2}\,dt \right| &\leq \int_0^\delta |h(t)| e^{-xt^2}\,dt \\
&< \int_0^\delta te^{-xt^2}\,dt \\
&< \int_0^\infty te^{-xt^2}\,dt \\
&= \frac{1}{2x}. \tag{5}
\end{align}
$$
Combining $(4)$ and $(5)$, $(3)$ becomes
$$
\int_0^\delta \frac{1}{1+t} e^{-xt^2}\,dt = \frac{1}{2} \sqrt{\frac{\pi}{x}} + O\!\left(\frac{1}{x}\right),
$$
and combining this with $(2)$ yields
$$
\int_0^\infty \frac{1}{1+t} e^{-xt^2}\,dt = \frac{1}{2} \sqrt{\frac{\pi}{x}} + O\!\left(\frac{1}{x}\right)
$$
for all $x > 1$.
Best Answer
Let's consider a bit more general case $$I(x,a)=\int_0^1\cos(xt^a)\tan tdt=\Re\int_0^1e^{ixt^a}\tan tdt\overset{xt^a=z}{=}\frac{1}{ax^{1/a}}\Re\int_0^xe^{iz}\tan\left(\frac{z}{x}\right)^\frac1az^{\frac1a-1}dz$$ Decomposing $\tan$ into the series $$=\frac{1}{ax^{1/a}}\Re\int_0^xe^{iz}\left(\big(\frac{z}{x}\big)^{1/a}+\frac{1}{3}\big(\frac{z}{x}\big)^{3/a}+...\right)z^{\frac 1a-1}dz\tag{1}$$ Let's consider the first term: $$I_1=\frac{1}{ax^{1/a}}\Re\int_0^xe^{iz}\Big(\frac{z}{x}\Big)^{1/a}z^{\frac 1a-1}dz=\frac{1}{ax^{2/a}}\Re\int_0^xe^{iz}z^{\frac2a-1}dz\tag{2}$$ Next, we integrate along a quarter-circle in the complex plane: $0\to x\to ix\to i0$, adding the arch of radius $x$. We do not have poles inside our contour; therefore, $$\oint_C=0\,\Rightarrow\,I_1+\frac{1}{ax^{2/a}}\Re\,e^{\frac{\pi i}{a}}\int_x^0e^{-z}z^{\frac2a-1}dz+I_{C_x}=0$$ where $I_{C_x}$ is the integral along the arch of the radius $x$. Extending integration of the second term to $\infty$ (and, therefore, dropping exponentially small corrections), $$I_1\sim\frac{\cos\frac{\pi}{a}\Gamma\left(\frac2a\right)}{ax^{2/a}}-I_{C_x}$$ To estimate $I_{C_x}$, we use the Jordan' lemma. Using $\sin\phi>\frac2\pi\phi$ for $\phi\in[0;\frac\pi2]$, $$|I_{C_x}|=\frac{1}{ax^{2/a}}\left|\int_0^{\pi/2}e^{ix\cos\phi-x\sin\phi}x^\frac2ae^{i\phi\frac2a}id\phi\right|<\frac{1}{a}\int_0^{\pi/2}e^{-\frac2\pi x\phi}d\phi<\frac{\pi}{2ax}$$ We see that $|I_{C_x}|\ll I_1$ at $x\to\infty$ and can be dropped , if $x^\frac2a\ll x$, or $\,\boxed{a>2}$.
In the same way we can evaluate next terms of $\tan$ decomposition and get, for example, $$I(x,a)=\int_0^1\cos(xt^a)\tan tdt\sim\frac{\cos\frac{\pi}{a}\Gamma\left(\frac2a\right)}{a\,x^{2/a}}+\frac{\cos\frac{2\pi}{a}\Gamma\left(\frac4a\right)}{3a\,x^{4/a}},\,\,a>4\tag{3}$$ For $a=3$ we are only allowed to keep the first term: $$\int^1_0 \cos(xt^3)\tan(t)dt\sim\frac{\cos\frac{\pi}{3}\Gamma\left(\frac23\right)}{3x^{2/3}}=\frac16\Gamma\big(\frac23\big)\,x^{-\frac23}\tag{4}$$
$\bf{Explicite \,form \,of\, full \,asymptotic}$
In his post @AxelT evaluated first oscillating asymptotic term, which gives very good approximation even at $x\sim1$. It would be also interesting to get full asymptotic.
Using integration in the complex plane, we got above $$I(x,a)=\frac{1}{ax^{1/a}}\Re\int_0^xe^{iz}\tan\left(\frac{z}{x}\right)^\frac1az^{\frac1a-1}dz=\frac{1}{ax^{1/a}}\Re \,e^\frac{\pi i}{2a}\int_0^xe^{-z}\tan\left(\frac{z^\frac1a}{x^\frac1a}e^\frac{\pi i}{2a}\right)z^{\frac1a-1}dz$$ $$-\,I_{C_x}\tag{5}$$ where $I_{C_x}$ is the integral along the arch of the radius $x$: $$I_{C_x}=\frac{1}{ax^{1/a}}\Re\int_0^{\pi/2}e^{ixe^{i\phi}}e^{\frac{i\phi}{a}}\tan\big(e^{\frac{i\phi}{a}}\big)x^\frac1ae^{i\phi}id\phi\overset{s=e^{i\phi}}{=}\frac1a\int_1^{e^{\frac{\pi i}{2}}}e^{ixs}\tan(s^\frac1a)s^\frac1ads\tag{6}$$ Integrating (6) several times by part and dropping exponentially small terms ($\sim e^{-\frac{\pi x}{2}}$) $$I_{C_x}\sim\frac{1}{a}\Re\left(-\frac{e^{ix}}{ix}f^{(0)}(s,a)|_{s=1}+\frac{e^{ix}}{(ix)^2}f^{(1)}(s,a)|_{s=1}-\frac{e^{ix}}{(ix)^3}f^{(2)}(s,a)|_{s=1}-...\right)$$ where we denoted $$f^{(k)}(s,a)|_{s=1}:=\frac{\partial^k}{\partial s^k}s^\frac1a\tan(s^\frac1a)\Big|_{s=1}$$ Rearranging $$-I_{C_x}\sim\frac{\sin x}{ax}\left(f^{(0)}(s,a)-\frac{f^{(2)}(s,a)}{x^2}+...\right)\bigg|_{s=1}+\frac{\cos x}{ax^2}\left(f^{(1)}(s,a)-\frac{f^{(3)}(s,a)}{x^2}+...\right)\bigg|_{s=1}$$ $$=\frac1a\sum_{n=1}^\infty\frac{\sin\big(\frac{\pi(n-1)}{2}+x\big)}{x^n}\left(s^\frac1a\tan s^\frac1a\right)^{(n-1)}\,\bigg|_{s=1}\tag{7}$$ Using the series of $\tan$ $$\tan t=\sum_{n=1}^\infty(-1)^{n-1}\frac{4^n(4^n-1)}{(2n)!}B_{2n}t^{2n-1}$$ and integrating term by term (extending integration to $\infty$) $$\frac{1}{ax^{1/a}}\Re \,e^\frac{\pi i}{2a}\int_0^xe^{-z}\tan\left(\frac{z^\frac1a}{x^\frac1a}e^\frac{\pi i}{2a}\right)z^{\frac1a-1}dz\sim\frac1a\sum_{n=1}^\infty\frac{(-1)^{n-1}4^n(4^n-1)B_{2n}}{(2n)!}\frac{\cos\frac{\pi n}{a}\Gamma\big(\frac{2n}{a}\big)}{x^\frac{2n}{a}}\tag{8}$$ Putting (7) and (8) into (5), we get the desired asymptotics: $$\color{blue}{I(x,a)\sim\frac1a\sum_{n=1}^\infty\frac{(-1)^{n-1}4^n(4^n-1)B_{2n}}{(2n)!}\frac{\cos\frac{\pi n}{a}\Gamma\big(\frac{2n}{a}\big)}{x^\frac{2n}{a}}}$$ $$\color{blue}{+\,\frac1a\sum_{n=1}^\infty\frac{\sin\big(\frac{\pi(n-1)}{2}+x\big)}{x^n}\left(s^\frac1a\tan s^\frac1a\right)^{(n-1)}\,\bigg|_{s=1}}$$ Here we have the only limitation $\color{blue}{a>0}$. If we want to take several first terms, the solution - the sequence of the terms we have to keep - depends, of course, on $a$. For example, for $a=3$ we get the series $$I(x,3)\sim\frac{\cos\frac\pi3\Gamma\big(\frac23\big)}{3x^\frac23}+\frac{\sin x\tan1}{3x}+\frac{\cos\frac{2\pi}{3}\Gamma\big(\frac43\big)}{9x^\frac43}+\frac{2\cos\pi\,\Gamma(2)}{45x^2}+\frac{\cos x}{9x^2}\left(\tan1+\frac{1}{\cos^21}\right)$$ $$+\,O\left(\frac1{x^\frac83}\right)$$